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ABSTRACT 

We perform a three-dimensional multi-probe analysis of the rich galaxy cluster A1689, one of the most pow¬ 
erful known lenses on the sky, by combining improved weak-lensing data from new wide-held BVRci 1 z' 
Subaru/Suprime-Cam observations with strong-lensing. X-ray, and Sunyaev-Zel’dovich effect (SZE) data sets. 

We reconstruct the projected matter distribution from a joint weak-lensing analysis of two-dimensional shear 
and azimuthally integrated magnification constraints, the combination of which allows us to break the mass- 
sheet degeneracy. The resulting mass distribution reveals elongation with an axis ratio of ~ 0.7 in projection, 
aligned well with the distributions of cluster galaxies and intracluster gas. When assuming a spherical halo, 
our full weak-lensing analysis yields a projected halo concentration of c|J/ 0c = 8.9 ± 1.1 (c 2 j r ’ ~ 11), consis¬ 
tent with and improved from earlier weak-lensing work. We find excellent consistency between independent 
weak and strong lensing in the region of overlap. In a parametric biaxial framework, we constrain the intrinsic 
structure and geometry of the matter and gas distributions, by combining weak/sbong lensing and X-ray/SZE 
data with minimal geometric assumptions. We show that the data favor a biaxial geometry with minor-major 
axis ratio 0.39 ± 0.15 and major axis closely aligned with the line of sight (22° ± 10°). We obtain a halo mass 
M2ooc = (1.2 ±0.2) x 10 15 Mq/i _1 and a halo concenbation C200C = 8.4 ± 1.3, which overlaps with the a) la- 
tail of the predicted distribution. The shape of the gas is rounder than the underlying matter but quite elongated 
with minor-major axis ratio 0.60 ± 0.14. The gas mass fraction within 0.9 Mpc is 10+2%, a typical value for 
high-mass clusters. The thermal gas pressure contributes to ~ 60% of the equilibrium pressure, indicating a 
significant level of non-thermal pressure support. When compared to Planck’ s hydrostatic mass estimate, our 
lensing measurements yield a spherical mass ratio of Mpi anc k/ Mgl = 0.70 ± 0.15 and 0.58 ± 0.10 with and 
without corrections for lensing projection effects, respectively. 

Keywords: cosmology: observations — dark matter — galaxies: clusters: individual (A1689) — gravitational 
lensing: weak — gravitational lensing: strong 


1. INTRODUCTION 

The evolution of the abundance of galaxy clusters with cos¬ 
mic epoch is sensitive to the amplitude and growth rate of pri¬ 
mordial density fluctuations as well as to the cosmic volume- 
redshift relation because massive clusters lie in the high-mass 
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exponential tail of the halo mas function (Haiman et al. 2001; 
Watson et al. 2014). Therefore, large cluster samples defined 
from cosmological surveys can provide an independent means 
of examining any viable cosmological model, including the 
current concordance A cold dark matter (ACDM) model de¬ 
fined in the framework of general relativity, complement¬ 
ing cosmic microwave background (CMB), large-scale galaxy 
clustering, and supernova observations. 

Clusters provide various probes of the role and nature of 
“dark matter” (DM) that dominates the material universe 
(Clowe et al. 2006), or modified gravity theories as an alter¬ 
native to DM (Rapetti et al. 2010), physics governing the final 
state of self-gravitating collisionless systems in an expanding 
universe (Navarro et al. 1996, 1997; Taylor & Navarro 2001; 
Hjorth & Williams 2010), and screening mechanisms in long- 
range modified models of gravity whereby general relativity 
is restored (Narikawa et al. 2013). 

Substantial progress has been made in recent years in con- 
sbucting statistical samples of clusters thanks to dedicated 
surveys (e.g., Planck Collaboration et al. 2014, 2015c; Bleem 
et al. 2015). Cluster samples are often defined by X-ray 
or Sunyaev-Zel’dovich effect (SZE) observables, so that the 
masses are indirectly inferred from scaling relations, which 
are often based on the assumption of hydrostatic equilibrium 
(HSE) and then statistically calibrated using weak lensing or 
internal dynamics using a subset of massive clusters at lower 
redshifts (Rines et al. 2013; Gruen et al. 2014). Since the 
level of mass bias from indirect observations assuming HSE 
is likely mass dependent (Sereno et al. 2014a) and sensitive 
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to calibration systematics of the instruments (Donahue et al. 
2014; Israel et al. 2015), a systematic effort is needed to en¬ 
able a self-consistent calibration of mass-observable relations 
using robust, direct cluster mass measurements (von der Lin¬ 
den et al. 2014; Umetsu et al. 2014; Merten et al. 2014; Ford 
et al. 2014; Jimeno et al. 2014; Hoekstra et al. 2015; Simet 
et al. 2015) and well-defined selection functions (e.g., Ben¬ 
itez et al. 2014). 

The great attraction of gravitational lensing in the cluster 
regime is its ability to map the mass distribution on an individ¬ 
ual cluster basis, independent of and free from assumptions 
about the physical and dynamical state of the cluster system 
(Miyazaki et al. 2007; Okabe & Umetsu 2008; Hamana et al. 
2009). Clusters act as efficient gravitational lenses, produc¬ 
ing various observable effects, including deflection, distor¬ 
tion, and magnification of the images of background sources 
(Bartelmann & Schneider 2001). In the weak regime, the lens¬ 
ing signals are approximately linearly related to the gravita¬ 
tional potential, so that one can determine the distribution of 
lensing matter at large scales in a model-independent manner 
(e.g., Umetsu et al. 1999, 2011b). In the strong regime, sev¬ 
eral sets of multiply-lensed images with known redshifts can 
be used to constrain the mass distribution in the cluster cores 
(e.g., Jauzac et al. 2014; Zitrin et al. 2014). 

A practical difficulty of obtaining precise mass measure¬ 
ments from cluster lensing, however, is significant scatter 
present in the projected lensing signals due to inherent varia¬ 
tions (at a fixed halo mass) in halo concentration, asphericity, 
orientation, and the presence of correlated large scale struc¬ 
ture (Rasia et al. 2012). The projection effects due to such 
intrinsic profile variations alone can produce a 1$ 20% un¬ 
certainty in lensing mass estimates for ~ 1O 15 M 0 clusters 
(Becker & Kravtsov 2011; Gruen et al. 2015). 

A possible way to overcome this problem is to simultane¬ 
ously determine the mass, concentration, shape, and orienta¬ 
tion of a given cluster by combining lensing data with inde¬ 
pendent probes or information about its line-of-sight elonga¬ 
tion (Sereno 2007; Corless et al. 2009; Limousin et al. 2013). 
Gravitational lensing probes the structure and morphology of 
the matter distribution in projection. X-ray observations con¬ 
strain the characteristic size and orientation of the intracluster 
medium (ICM) in the sky plane. The elongation of the ICM 
along the line of sight can be constrained from the combina¬ 
tion of X-ray and thermal SZE observations (De Filippis et al. 
2005; Sereno et al. 2012). Recently, Sereno et al. (2013) de¬ 
veloped a parametric triaxial framework to combine and cou¬ 
ple independent morphological constraints from lensing and 
X-ray/SZE data, using minimal geometric assumptions about 
the matter and gas distributions but without assuming HSE. 

The first critical step in a three-dimensional (3D) cluster 
analysis is an unbiased, direct recovery of the projected clus¬ 
ter mass distribution from weak lensing. A fundamental lim¬ 
itation of measuring shear only is the mass-sheet degeneracy 
(Schneider & Seitz 1995). This degeneracy can be broken by 
using the complementary combination of shear and magnifi¬ 
cation (Schneider et al. 2000; Umetsu et al. 2011b; Umetsu 
2013). Umetsu et al. (2011b) have shown that the magnifi¬ 
cation effect can significantly enhance the accuracy and pre¬ 
cision of lensing-derived cluster mass profiles when added to 
weak-lensing shear measurements. 

Our aim in this paper is to develop and apply a comprehen¬ 
sive set of techniques and methods for 3D analysis of galaxy 
clusters based on the multi-probe framework of Sereno et al. 


Table 1 

Properties of the galaxy cluster A1689 


Parameter 

Value 

ID. 

A1689 

Optical center position (J2000.0) 

R.A. 

13:11:29.52 

Decl. 

-01:20:27.59 

X-ray center position (J2000.0) 

R.A. 

13:11.29.50 

Decl. 

-01:20:29.92 

SZE center position (J2000.0) 

R.A. 

13:11.29.57 

Decl. 

-01:20:29.87 

Redshift. 

0.183 

X-ray temperature (keV). 

10.4 

Einstein radius (") . 

47.0 ± 1.2 atz 3 = 2 


References. — [1] Andersson & Madejski (2004); [2] 
Limousin et al. (2007); [3] Kawaharada et al. (2010); [4] 
Coe et al. (2010). 


Note. — The optical cluster center is defined as the 
center of the BCG from Ref. [2], Units of right ascension 
are hours, minutes, and seconds, and units of declination 
are degrees, arcminutes, and arcseconds. The X-ray prop¬ 
erties were taken from Ref. [3]. The X-ray center is de¬ 
fined as the X-ray emission centroid derived from XMM- 
Newton observations. See also Ref. [1], The SZE cen¬ 
ter is determined from the joint analysis of interferometric 
BIMA/OVRO/SZA observations described in Section 7.2. 

The BCG is located within 2. 3" (— 5kpc h ~ 1 ) of the X- 
ray center. The X-ray and SZE centroid positions agree to 
within 1". The Einstein radius is constrained by detailed 
strong lens modeling by Ref. [4], 

(2013). For this aim, we first generalize the one-dimensional 
(ID) weak-lensing inversion method of Umetsu et al. (201 lb) 
to a two-dimensional (2D) description of the mass distribution 
without assuming particular functional forms, i.e., in & free¬ 
form fashion. In this approach, we combine the spatial shear 
pattern with azimuthally averaged magnification information, 
imposing integrated constraints on the mass distribution. 

Taking advantage of new BVRci'z' imaging obtained with 
Suprime-Cam on the 8.3 m Subaru Telescope, we perform 
a new weak-lensing analysis of the rich cluster A1689 at 
z = 0.183 and then apply our methods to weak-lensing, 
strong-lensing. X-ray, and SZE data sets we have obtained 
for the cluster. The cluster is among the best studied clus¬ 
ters (Tyson & Fischer 1995; Taylor et al. 1998; Andersson & 
Madejski 2004; Broadhurst et al. 2005b; Halkola et al. 2006; 
Limousin et al. 2007; Umetsu & Broadhurst 2008; Peng et al. 
2009; Kawaharada et al. 2010; Coe et al. 2010; Sereno et al. 
2012; Nieuwenhuizen & Morandi 2013; Sereno et al. 2013) 
and one of the most powerful known lenses on the sky, char¬ 
acterized by a large Einstein radius of (?Ein = 47.0" ±1.2" for 
a fiducial source at t; s = 2 (see Table 1; Coe et al. 2010); this 
indicates a high degree of mass concentration in projection 
(Broadhurst & Barkana 2008). To date, 61 candidate systems 
of 165 multiply-lensed images have been identified (Broad¬ 
hurst et al. 2005b; Coe et al. 2010; Diego et al. 2015) from Ad¬ 
vanced Camera for Surveys (ACS) observations with the Hub¬ 
ble Space Telescope (HST). Despite significant efforts, the de¬ 
gree of concentration inferred from different lensing analyses 
is somewhat controversial (see Coe et al. 2010; Sereno et al. 
2013), and it is still unclear if and to what degree this cluster 
is over-concentrated. 

The paper is organized as follows. After summarizing the 
basic theory of cluster weak lensing, we present in Section 
2 the formalism that we use for our weak-lensing analysis. 
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In Section 3, we describe our Subaru observations and data 
processing. In Section 4, we present our Subaru weak-lensing 
analysis. Section 5 presents our HST strong-lensing analysis. 
In Section 6 we outline the triaxial modeling and describe the 
statistical framework for the 3D cluster analysis. In Section 7 
we present the multi-probe analysis of lensing and X-ray/SZE 
data. In Section 8 we discuss the results and their implications 
for the intrinsic properties of A1689. Finally, a summary of 
our work is given in Section 9. 

Throughout this paper, we use the AB magnitude system 
and adopt a concordance ACDM cosmology with f l m = 
0.3, 17a = 0.7, and h = 0.7h^ o = 0.7 where H 0 = 
h x lOOkms -1 Mpc -1 . In this cosmology, V corresponds 
to 129 kpc /i -1 ~ 185 kpc /iT ,, 1 for this cluster. The refer¬ 
ence sky position is the center of the brightest cluster galaxy 
(BCG): R.A.(J2000.0) = 13 : 11 : 29.52, Decl.(J2000.0) = 
—01 : 20 : 27.59 (Table 1). We use the standard notation 
/'a to denote the spherical overdensity radius within which 
the mean interior density is A times the critical density p c of 
the universe at the cluster redshift. For its ellipsoidal counter¬ 
part R/\, see Section 6.1. All quoted errors are 68.3% (lcr) 
confidence limits (CF) unless otherwise stated. 


2. WEAK-LENSING METHODOLOGY 

2.1. Weak Lensing Basics 

In the cluster regime, the lensing convergence, k = E/E c , 
is the projected mass density E(0) in units of the critical 
surface density for lensing, E c = (c 2 D s )/{M:GD\D\ S ) = 
c 2 /(AirGDiP) with D\, D s , and D\ s the lens, source, and 
lens-source angular diameter distances, respectively; f3(z) = 
Di s (z)/D s (z) represents the geometric lensing strength for a 
source at redshift z, where /3(z) = 0 for 2 < z\. 

The gravitational shear 7 = 71 + 172 can be directly ob¬ 
served from ellipticities of background galaxies in the weak 
regime, I. The shear and convergence are related by 

7(0) = fM D(e-*)«&) (i) 

with 79(0) = ( 6 \ — 6 \ — 216162 )/{it\0\^) (Kaiser & Squires 
1993). The observable quantity for quadrupole weak lensing 
in general is not 7 but the complex reduced shear. 


9 { 0 ) 


7(A) 

1 — k(6)' 


( 2 ) 


The g field is invariant under k(0) —> Xk(6) -+■ 1 — A and 
7 ( 0 ) -A X'y(O) with an arbitrary constant A / 0, known as 
the mass-sheet degeneracy (Schneider & Seitz 1995). This de¬ 
generacy can be broken, for example, by measuring the mag¬ 
nification /t( 0 ) in the subcritical regime, 

/i(0)= [1-«(0)P-| 7 (0)I 2 “ A(0)’ P) 

which transforms as p(8) -A A 2 p(Q). 

Fet us consider a population of source galaxies described by 
their redshift distribution function, N(z). In general, we ap¬ 
ply different size, magnitude, and color cuts in background se¬ 
lection for measuring shear and magnification, which results 
in different N(z). In contrast to the former effect, the latter 
does not require source galaxies to be spatially resolved, but it 
requires a stringent flux limit against incompleteness effects. 


The mean lensing depth for a given population (X = g. //) is 


(p)x = 

■ /»oo 

/ dzw(z)N x {z)j3(z) 

■ poo 

/ dzw{z)N x{z) 


Jo 

Jo 

( 


where w(z) is a weight factor (see Section 3.3). 

We introduce the relative lensing strength of a given source 
population relative to a fiducial source in the far background 
as (IF) a' = {P)x/Poo (Bartelmann & Schneider 2001) with 
Ax, = /3(z -A 00 ; Zi). The associated critical density is 
X Ci oo(-Zi) = c 2 /(AnGDipoo). Hereafter, we use the far- 
background fields Hoc ( 0 ) and 700 ( 0 ) to describe the projected 
cluster mass distribution. 


2.2. Discretized Mass Distribution 
We discretize the convergence field n oo (0) = E“J o X(0) 
into a regular grid of pixels and approximate (0 ) by a lin¬ 
ear combination of basis functions B(0 — Q') as 

N pix 

Koo (0) = E-J o ^H(0-0 ro )E m , (5) 

m—1 

where our model (signal) s = { E ,„ is a vector of param¬ 
eters containing mass coefficients. To avoid the loss of infor¬ 
mation due to oversmoothing, we take the basis function to be 
the Dirac delta function 77(0 —0 m ) = (A0) 2 <5p(0 — 0 m ) with 
A 9 a constant spacing, so that s represents the cell-averaged 
projected mass density. The 7 ^ (0) field can be expressed as 

tV pix 

7oo(0) = E V (° - (6) 

m—1 

with V = D ig) B an effective kernel (Equation (1)). Hence, 
both k x and 7 ^ can be written as linear combinations of s. 

Because of the choice of the basis function, an unbiased ex¬ 
traction of mass coefficients {Eml^jEj (or certain linear com¬ 
binations of E m ) can be done by performing a spatial integral 
of Equation (5) over a certain area. In practical applications, 
such operations include smoothing (Figure 1), azimuthal av¬ 
eraging for a mass profile reconstruction (Section 5.3), and 
profile fitting with smooth functions (Section 7). 


2.3. Weak-lensing Observables 
2.3.1. Reduced Shear 


The quadrupole image distortion due to lensing is described 
by the reduced shear, g = gi +ig 2 - We calculate the weighted 
average g m = g(0 m ) of individual shear estimates on a regu¬ 
lar cartesian grid (m = 1 , 2 ,..., 7V p i x ) as 


. k 

(7) 

where *S'(0(fe), 0m.) is a spatial window function, g^) is an 
estimate of g{6) for the A;th object at 0 (fc), and w^k) is its sta¬ 
tistical weight given by w^ k) = 1 /(<r 2 (fc) + a 2 ) with cr 2 (fc) the 
error variance of g^ k ) and a 2 the softening constant variance. 

We choose a g = 0.4, a typical value of the mean rms \J~^ 
found in Subaru observations (e.g., Umetsu et al. 2009). 


9m — 


5 


. k 
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The source-averaged theoretical expectation for the estima¬ 
tor (7) is approximated by (see Appendix A.l) 


9 ( 9 m ) 


(W) g 7oo (e m ) 

1- fwAW) g ^(9my 


( 8 ) 


where (W) g is the source-averaged relative lensing strength 
(Section 2.1), and fw, g = (W 2 ) g /(W ) 2 is a dimensionless 
quantity of the order unity. The variance a 2 m = o 2 (6 m ) for 
g m = gi,m + ig 2 ,m is expressed as 


2011; Ford et al. 2012, 2014). The faint blue population ly¬ 
ing at z ~ 2 (e.g., Lilly et al. 2007; Medezinski et al. 2010, 
2013) tends to have a steep intrinsic slope close to the lensing- 
invariant one, a = 1. 

The covariance matrix of Ng(9) includes both sample co- 
variance and Poisson variance (Hu & Kravtsov 2003): 

Cov[N^e m ),N^e n )] = (C N ) mn = (N») 2 u mn +5 mn N^0 m ), 

( 13 ) 

where u; mn is the cell-averaged angular correlation function 


a 


2 

g,m 


& (®(fc)> @rn)'W(k) (J g(k) 

_ k 


^ ' '5 (9(k) i 9m)W(k) 

_ k 


( 9 ) 

In this work, we adopt the top-hat window of radius 9f 
(Merten et al. 2009), S(9,0') = H(9 t - \0 - 0'\), with H(x) 
the Heaviside function defined such that H(x) = 1 if x > 0 
and H (x) = 0 otherwise. The covariance matrix for g m is 


-2 


Gov(g a ^m, gg,n) — 8 a {3 (Gg) rnn — ~ tT grn Ct gn ^H ( 1 9 m $n|); 

( 10 ) 


where £,n(x\ 9{) is the autocorrelation of a pillbox of radius 
9{ (White et al. 1999; Park et al. 2003), given by 


uJmn = U 2 ed 2 e' s m (e)s n (0')u(0-0') ( 14 ) 

* *cell J 

with u>(0) the angular two-point correlation function of the 
source galaxies, S m (9) the boxcar window function of the 
nvth cell, and H ce n = fd 2 9 S m (0). For deep lensing obser¬ 
vations, the angular correlation length of background galax¬ 
ies can be small (e.g., Connolly et al. 1998) compared to the 
typical resolution ~ 1' of reconstructed mass maps. There¬ 
fore, the correlation between different cells can be generally 
ignored, whereas the unresolved correlation on small angular 
scales accounts for increase of the variance of N fl (d) (Van 
Waerbeke et al. 2000). We thus approximate Cy by 

(Gn) rnn ~ [<<W M (0 ro )> +iV M (0 m )] S mn , (15) 



(ID 

for |a?| < 26f and £h(x) = 0 for \x\ > 29{. 


with (8N 2 (0 m )) the variance of the mth counts. 

To enhance the signal-to-noise ratio, we azimuthally aver¬ 
age Ng,(d) in contiguous, concentric annuli and calculate the 
surface number density of background galaxies as 

a function of clustercentric radius: 


2.3.2. Magnification Bias 

Deep multi-band photometry allows us to explore the faint 
end of the luminosity function of red quiescent galaxies at 
z ~ 1 (Ilbert et al. 2010), for which the effect of magnification 
bias is dominated by the geometric area distortion and thus not 
sensitive to the exact form of the source luminosity function. 
In this work, we perform magnification measurements using 
a flux-limited sample of red background galaxies. 

If the magnitude shift 8m = 2.5 log 10 ji due to magnifica¬ 
tion is small compared to that on which the logarithmic slope 
of the luminosity function varies, their number counts can 
be locally approximated by a power law at the limiting flux 
(Broadhurst et al. 1995). The expectation value for the source 
counts N /J ,(6 rn ) on a grid of equal-area cells (m = 1,2,...) is 
modified by lensing magnification as (see Appendix A. 2) 

E[N^(O m )\ =NgA 1 ~ a (6 m ), 

A(d) = [1 - (W), Koo (e)} 2 (W)l\ loc (9)\ 2 , 

where N M is the unlensed mean source counts per cell, a 
is the unlensed count slope evaluated at the flux limit F, 
a = — dlogA r /i (> F)/dlog.F, 14 and (W)^ is the source- 
averaged relative lensing strength (Section 2.1). 

The net magnification effect on the source counts vanishes 
when a = 1. In the regime where a < 1, the bias is domi¬ 
nated by the expansion of the sky area, producing a net count 
depletion. For a population with a > 1, the bias is positive, 
and a net density enhancement results (e.g., Hildebrandt et al. 

14 In the weak-lensing literature, s = dlog 10 N(< m)/dm = 0.4a in 
terms of the limiting magnitude m is often used instead of a (e.g., Umetsu 
et al. 2011b, 2014; Medezinski et al. 2013). 


n ^= <r-Y l ' p imN lt (0 ra ) ( 16 ) 

with Vim = (Em Ami) -1 Ami the radial projection matrix 
normalized as Yhm Ptm = 1. Here A m , represents the frac¬ 
tion of the area of the mth cell lying within the ith annular bin 
(0 < A mi < 1), and 7?i(> 1) is the mask correction factor 

for the tth annular bin, ru = E m (! - /m)A mi ] _1 A mi, 
with /„, the fraction of the mask area in the mth cell, due to 
bad pixels, saturated objects, foreground and cluster member 
galaxies (see Section 3.2 of Umetsu et al. 2014). 

The theoretical expectation for the estimator (16) is 

h^ i = n ll Y,' p imA 1 - a (9m) ( 17 ) 


with rig, = Ng/Q, ce \\. The bin-to-bin covariance matrix for 
the estimator (16) is obtained as 


Co v(ng t i,ngj) — (Cg) — 2 2 Vj m Vjn (Cjv) 

“cell 


( 18 ) 


Note that since Cn is diagonal, Gj, is also diagonal: 

(Cg)ij = Gggdij- ( 19 ) 


2.4. Mass Reconstruction 

Given a model m and observed (fixed) data <7, the posterior 
probability P(m\d) is proportional to the product of the like¬ 
lihood C(m) = P(d\m) and the prior probability P(m). In 
our 2D inversion problem, m is a vector containing the sig¬ 
nal parameters s (Section 2.2) and calibration parameters c 
(Section 2.4.3), m = ( s,c ). 
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The total likelihood function £ for combined weak-lensing 
data d is given as a product of the two separate likelihoods, 
£ = CgC-n , where £,, and £ p are the likelihood functions 
for shear and magnification, respectively. We assume that the 
errors on the data follow a Gaussian distribution, so that £ ex 
exp(—X 2 /2), with \ 2 the standard misfit statistic. 

2.4.1. Shear Log-likelihood Function 

The log-likelihood function l g = — In £ g for 2D shear data 
can be written in the general form (ignoring constant terms) 
as (Oguri et al. 2010; Umetsu et al. 2012) 

1 AW 2 

lg = 2 5Z 0^'g)mn [9ot9a,n{'m)\ 

m,n— 1 a—1 

( 20 ) 

where g a ,m(tn) is the theoretical expectation for g a ,m = 
g a (6 m ), and (VV g ) rra „ is the shear weight matrix. 


Table 2 

Subaru/Suprime-Cam data 


Filter 

Exposure time 3 
(ks) 

Seeing b 
(arc sec) 

mum 0 
(AB mag) 

B 

2.40 

0.91 

27.1 

V 

4.08 

0.84 

27.0 

Rc 

6.42 

0.70 (0.60) 

27.0 

i' 

4.08 

0.84 

26.4 

z' 

8.02 

0.81 

26.2 


Note. — The Rq band is used as the filter to mea¬ 
sure object shapes for the weak-lensing analysis, where 
we separately stack data from different epochs. The 
J?C"band seeing in parentheses is the average of val¬ 
ues derived from separate stacks. 
a Total exposure time. 

b Seeing FWHM in the full stack of images. 
c Limiting magnitude for a 3cr detection within a 2" 
aperture. 


(Wg) mn = M mM n (C^) mn , ( 21 ) 

with (C~ 1 ) mn the inverse covariance matrix for the 2D shear 
data and M m a mask weight, defined such that M m = 0 if the 
mth cell is masked out and M m = 1 otherwise. 

2.4.2. Magnification Log-likelihood Function 

Similarly, the log-likelihood function for magnification- 
bias data l p = — In £ p can be written as 

1 AW 

If- = 2 y DViJ ~ (22) 

i=l 

where n p ti (m) is the theoretical prediction for the observed 
counts n p> i (see Equations (16) and (17)), and is the 

magnification weight matrix, 

= ( c £)„ = ®) 
u ii,i 

(Equations (18) and (19)). We use Monte Carlo integration 
to calculate the radial projection matrix Vi m (Equation 16) of 
size lVbi n x TVpix, which is needed to predict 
for a given m = ( s , c). 

The / ( , function imposes a set of azimuthally integrated con¬ 
straints on the underlying projected mass distribution. Since 
magnification is locally related to k, this will essentially pro¬ 
vide the otherwise unconstrained normalization of E (6) over 
a set of concentric rings where count measurements are avail¬ 
able. We note that no assumption is made of azimuthal sym¬ 
metry or isotropy of the 2D mass distribution E (0). 

2.4.3. Calibration Parameters 

We account for the calibration uncertainty in the observa¬ 
tional nuisance parameters, 

c={(W) g J w , g ,(W)^n ll ,a). (24) 

To do this, we include in our analysis Gaussian priors on c 
given by means of quadratic penalty terms with mean values 
and errors directly estimated from data. 

2.4.4. Best-fit Solution and Covariance Matrix 

The log posterior F(m) = — lnP(m|d) is expressed 
as a linear sum of the log-likelihood and prior terms. The 


maximum-likelihood (ML) solution, rh, is obtained by mini¬ 
mizing F{m) with respect to m. In our implementation we 
use the conjugate-gradient method (Press et al. 1992) to find 
the solution. Here we employ an analytic expression for the 
gradient function V F(rn) obtained in the nonlinear, subcrit- 
ical regime. To be able to quantify the errors on the recon¬ 
struction, we evaluate the Fisher matrix atm = rh, as 

d 2 F(m) 

dmpdnipi 

where the angular brackets represent an ensemble average, 
and the indices (p. p r ) run over all model parameters m = 
(s, c). We estimate the error covariance matrix as 

Cov(m p , m p ') = Cp P ’ = . (26) 

3. SUBARU OBSERVATIONS 

Here we present a description of our data analysis of A1689 
based on deep Subaru BVRqi' z 1 images. In this study, we 
analyze the data using the same methods and procedures as in 
Umetsu et al. (2014), who performed a weak-lensing analy¬ 
sis of 20 high-mass clusters selected from the CLASH survey 
(Postman et al. 2012). For details of our reduction and analy¬ 
sis pipelines, we refer to Section 4 of Umetsu et al. (2014). 

3.1. Data and Photometry 

We analyze deep BVRci'z' images of A1689 observed 
with the wide-held camera Suprime-Cam (34' x 27'; Miyazaki 
et al. 2002) at the prime focus of the 8.3 m Subaru Telescope. 
We combine both existing archival data taken from SMOKA 15 
with observations acquired by the team on the nights of 2010 
March 17-18 (S10A-019). The observation details of A1689 
are summarized in Table 2. Figure 1 shows a BVRci'z' com¬ 
posite color image of the cluster held, produced using the pub¬ 
licly available TRILOGY software (Coe et al. 2012). The im¬ 
age is overlaid by mass contours determined from our weak- 
lensing analysis (see Section 4.2). 

Our imaging reduction pipeline derives from Nonino et al. 
(2009) and has been optimized separately for accurate pho¬ 
tometry and shape measurements. For multi-band photome¬ 
try, standard reduction steps include bias subtraction, super¬ 
flat-held correction, and point-spread-function (PSF, here- 

15 



(25) 


http://smoka.nao.ac.jp 
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Figure 1. Subaru BVRci' z' composite color image centered on the galaxy cluster A1689 (z = 0.183), overlaid with mass contours from our joint shear- 
and-magnification weak-lensing analysis of Subaru data. The image is 30 / X 30 ; in size. The mass map is smoothed with a Gaussian of FWHM = 1.5'. 
The horizontal bar represents 1 Mpc h ~ 1 at the cluster redshift. The lowest contour level and the contour interval are Ak = 0.06. The cyan contours show the 
smoothed projected distribution of cluster red-sequence galaxies. 

North is up and east is to the left. 


after) matching between exposures in the same band. An ac¬ 
curate astrometric solution is derived with the SCAMP soft¬ 
ware (Bertin 2006), using the the Sloan Digital Sky Survey 
(SDSS, Adelman-McCarthy et al. 2008) as an external refer¬ 
ence catalog. 16 The Swarp software (Bertin et al. 2002) is 
used to stack individual exposures on a common World Coor¬ 
dinate System (WCS) grid with pixel scale of 0.2". 

The photometric zero-points for the co-added images were 
derived using //5T/ACS magnitudes of cluster elliptical-type 
galaxies. These zero points were further refined by fit- 

16 This research has made use of the VizieR catalog access tool, CDS, 
Strasbourg, France. 


ting SED (spectral energy distribution) templates with the 
BPZ code (Bayesian photometric redshift estimation; Benitez 
2000; Benitez et al. 2004) to 1445 galaxies having spectro¬ 
scopic redshifts. 17 This leads to a final photometric accu¬ 
racy of ~ 0.01 mag in all passbands. The magnitudes were 
corrected for Galactic extinction according to Schlegel et al. 
(1998). The multi-band photometry was measured using SEx- 
tractor (Bertin & Arnouts 1996) in dual-image mode on PSF- 
matched images created by ColorPro (Coe et al. 2006). 

17 The data used here are part of an extensive multi-object spectroscopy 
survey carried out with the VIMOS spectrograph on the VLT (Czoske 2004). 
For details, see Lemze et al. (2009). 
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3.2. Shape Measurement 

We use our shear analysis pipeline based on the IMCAT 
package (Kaiser et al. 1995, KSB) incorporating improve¬ 
ments developed by Umetsu et al. (2010). On the basis of 
simulated Subaru/Suprime-Cam images (Oguri et al. 2012; 
Massey et al. 2007), Umetsu et al. (2010) showed that the 
lensing signal can be recovered with \m\ ~ 5% of the multi¬ 
plicative shear calibration bias (as defined by Heymans et al. 
2006; Massey et al. 2007), and c ~ 10 -3 of the residual shear 
offset, which is about one order of magnitude smaller than 
the typical shear signal in cluster outskirts. Accordingly, we 
include for each galaxy a shear calibration factor of 1/0.95 
(g —>• <7/0.95) to account for residual calibration. 

In this work, we perform weak-lensing shape analysis us¬ 
ing the same procedures adopted in the CLASH weak-lensing 
analysis of Umetsu et al. (2014). Here, we only highlight key 
aspects of our analysis pipeline: 

• Object detection. Objects are detected using the IM¬ 
CAT peak finder, hfindpeaks , using a set of Gaussian 
kernels of varying sizes. This algorithm produces ob¬ 
ject parameters such as the peak position, the best- 
matched Gaussian scale length, r g , and an estimate of 
the significance of the peak detection, v. 

• Crowding effects. Objects having any detectable neigh¬ 
bors within 3 r g are identified. All such close pairs of 
objects are rejected to avoid possible shape measure¬ 
ment errors due to crowding. The detection threshold 
is set to v = 7 for close-pair identification. After this 
close-pair rejection, objects with low detection signifi¬ 
cance v < 10 are excluded from our analysis. 

• Shear calibration. We calibrate KSB’s isotropic correc¬ 
tion factor P g as a function of object size (r g ) and mag¬ 
nitude, using galaxies detected with high significance 
v > 30 (Umetsu et al. 2010). This is to minimize the 
inherent shear calibration bias in the presence of noise. 
We correct for the isotropic smearing effect caused by 
seeing as well as by the window function used in the 
shape estimate as g a = e a /P g with e a the anisotropy- 
corrected object ellipticity. 

To measure the shapes of background galaxies, we use the 
Ii c band data, which have the best image quality in our data 
sets (Table 2). Two separate co-added Ac-hand images are 
created, one from 2009 (observed by Matsuda et al.) and 
another from 2010 (observed by Umetsu et al.). We sep¬ 
arately stack data obtained at different epochs. We do not 
smear individual exposures before stacking, so as not to de¬ 
grade the weak-lensing signal. After PSF anisotropy correc¬ 
tion, the mean residual stellar ellipticity is consistent with 
zero, and the rms residual stellar ellipticity in each stack is 
er(<5e*) ~ 2.5 x 10 -3 per component. A shape catalog is 
created for each epoch separately. These subcatalogs are then 
combined by properly weighting and stacking the calibrated 
shear estimates for galaxies in the overlapping region (see 
Section 4.3 of Umetsu et al. 2014). 

3.3. Background Galaxy Selection 

A careful background selection is critical for a cluster 
weak-lensing analysis, so that unlensed objects do not dilute 
the true lensing signal of the background (Medezinski et al. 
2007; Umetsu & Broadhurst 2008; Okabe et al. 2013; Hwang 


Table 3 

Background Galaxy Samples for Weak-lensing Shape 
Measurements 


Sample 


n s a 

(arcmin ) 

2 eff b 

(D ls /D s ) 

fw 

Red 

12674 

12.0 

1.10 

0.79 ± 0.04 

1.00 

Blue 

9238 

8.7 

1.62 

0.84 ± 0.04 

1.01 

Blue+red 

21912 

20.7 

1.22 

0.80 ± 0.04 

1.01 


a Mean surface number density of source background galaxies. 
b Effective source redshift corresponding to the mean lensing depth 

(/3) = (D ls /D s ), defined as /3(z eB ) = (/3). 


2.5 



Rc~z' 


Figure 2. “Blue” and “red” background galaxy samples selected for the 
weak-lensing analysis (lower-left blue and lower-right red regions, respec¬ 
tively) on the basis of Subaru BRqz' color-color-magnitude selection. All 
galaxies (cyan) are shown in the diagram. At small clustercentric radius 
(< 4/), an overdensity of cluster galaxies is identified as our “green” sample 
(green), comprising mostly the red sequence of cluster ellipticals and some 
blue trail of later-type cluster members. The background samples are well 
isolated from the green region and satisfy other criteria as discussed in Sec¬ 
tion 3.3. The black dots represent a dynamically-selected spectroscopic sam¬ 
ple of 377 cluster galaxies found within a projected distance of 12' (~ r 200 c) 
from the cluster center. Our background selection successfully excludes all 
except 2 spectroscopically confirmed cluster members (see Section 3.3). 

et al. 2014). In particular, dilution due to contamination by 
cluster members can lead to a substantial underestimation of 
the true signal at small cluster radii, r r25ooc (Medezinski 
et al. 2010; Okabe et al. 2010). The relative importance of the 
dilution effect indicates that, the impact of background purity 

__^ / 2 

and depth is more important than that of shot noise (oc n g ). 

We use the color-color (CC) selection method of Medezin¬ 
ski et al. (2010) to define uncontaminated samples of back¬ 
ground galaxies from which to measure the shear and magni¬ 
fication effects. Here we refer the reader to Medezinski et al. 
(2010) for further details. Our multi-color approach and its 
variants have been successfully applied to a large number of 
clusters (Medezinski et al. 2010, 2011, 2013; Umetsu et al. 
2010, 2011b, 2012, 2014; Coe et al. 2012; Oguri et al. 2012; 
Covone et al. 2014; Sereno et al. 2014b). 

We use the Subam BRcz' photometry, which spans the full 
optical wavelength range, to perform CC selection of back¬ 
ground samples. In Figure 2, we show the B — Rq versus 
Rc — z' distribution of all galaxies to our limiting magnitudes 
(cyan). We select two distinct populations that encompass the 
red and blue branches of background galaxies in CC space, 
each with typical redshift distributions peaked around z ~ 1 
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Table 4 

Background Galaxy Samples for Magnification-bias Measurements 


Sample 

/ a 
^cut 

(AB mag) 


(arcmin -2 ) 

Q C 

<z> d 

Zeff e 

(■ D U /D S ) 

Red 

25.6 

26136 

19.0 ±0.5 

0.39 ± 0.08 

1.13 

1.05 

0.73 ± 0.04 

Blue 

25.6 

12143 

8.8 ±0.3 

0.82 ± 0.12 

1.81 

1.39 

0.82 ± 0.04 


a Fainter magnitude cut of the background sample. Apparent magnitude cuts are applied in the 
reddest CC-selection band available (z') to avoid incompleteness near the detection limit. 
b Coverage- and mask-corrected normalization of unlensed background source counts. 
c Logarithmic slope of the unlensed source counts a = 2.5 Tdlogin Nn(< z')/dz'] , . 

d Mean photometric redshift of the sample obtained with the BPZ code, defined similarly to 
Equation (4). 

e Effective source redshift corresponding to the mean lensing depth (/3) = (D\ s /D s ), defined 

as/3(2 e ff) = (/3>- 
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Figure 3. Azimuthally averaged radial profiles of the tangential lens distor¬ 
tion g. |_ (upper panel) and the 45° rotated (x) component g x (lower panel) 
for our red (triangles), blue (circles), and blue+red (squares) galaxy samples 
derived from Subaru multi-color photometry (Table 3). 

and ~ 2, respectively (see Figures 5 and 6 of Medezinski et al. 
2011; Lilly et al. 2007). The color boundaries of our “blue” 
and “red” background samples are shown in Figure 2. 

As a cross-check we calculate the tangential ( g + ) and cross 
( g x ) reduced-shear components in clustercentric radial bins, 
which we show in Figure 3. In the absence of higher-order 
effects, weak lensing produces only curl-free tangential dis¬ 
tortions, g + . The presence of x modes can thus be used to 
check for systematic errors. Using the weak-lensing-matched 
blue and red samples, we find a consistent, rising distortion 
signal all the way to the cluster center. For all cases, the 
x-component is consistent with a null signal detection well 
within 2cr at all radii. 

For the number counts to measure magnification, we de¬ 
fine flux-limited photometry samples of background galaxies. 
Here we limit the data to z l = 25.6 mag in the reddest band 
(Table 4), corresponding to the 5 a limiting magnitude within 
2" diameter aperture. We plot in Figure 4 the coverage- and 
mask-corrected surface number density as a function of clus¬ 
tercentric radius, for the blue and red samples. No cluster¬ 
ing is observed toward the center, demonstrating that there is 
no detectable contamination by cluster members in the back¬ 
ground samples. The red sample reveals a systematic decrease 
in their counts toward the cluster center, caused by magnfica- 
tion of the sky area (Section 3.3). The faint blue counts, on 


Figure 4. Coverage- and mask-corrected surface number density profiles 
of Subaru BRqz '-selected galaxy samples (Table 4). The results are shown 
for our red (triangles) and blue (circles) background samples. The error bars 
include contributions from Poisson counting uncertainties and contamination 
due to intrinsic clustering of each source population. For the red sample, a 
systematic radial depletion of the source counts is seen toward the cluster 
center owing to magnification of the sky area, while the faint blue counts are 
nearly constant with the distance from the cluster center. See also Figure 5. 

the other hand, are nearly constant with cluster radius, as ex¬ 
pected by their steep count slope (Table 4). A more quantita¬ 
tive magnification analysis is given in Section 4.1. 

For validation purposes, we compare in Figure 2 our back¬ 
ground samples with a dynamically-selected spectroscopic 
sample of 377 cluster galaxies (black) found within a pro¬ 
jected distance of 12' (~ r2ooc) from the cluster center. We 
find that our background selection procedure successfully ex¬ 
cludes all except 2 spectroscopically confirmed cluster mem¬ 
bers (see also Coe et al. 2012; Umetsu et al. 2012), corre¬ 
sponding to a negligible contamination fraction of ~ 0.5%. 
We note that, in the blue background region, there are 4 clus¬ 
ter members, of which two are excluded by the magnitude 
cuts used to reject bright foreground/cluster galaxies. 

We estimate the mean depths ((/?}, (/3 2 )) of the background 
samples (Tables 3 and 4), which are necessary when convert¬ 
ing the observed lensing signal into physical mass units. For 
this, we follow the prescription outlined in Section 4.4 of 
Umetsu et al. (2014). We utilize BPZ to measure photo-xs 
using our PSF-corrected Subam BVRci'z 1 photometry. Fol¬ 
lowing Umetsu et al. (2012), we employ BPZ’s ODDS param- 
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eter as the weight factor w(z) in Equation (4). The resulting 
depth estimates are summarized in Tables 3 and 4. 

4. SUBARU WEAK-LENSING ANALYSIS 

We use our z'-band limited sample of red background 
galaxies (Table 4) for magnifciation measurements and a full 
composite sample of blue+red galaxies (Table 3) for shear 
measurements. In Section 4.1, we perform a ID weak-lensing 
analysis of A1689 to derive azimuthally averaged lensing pro¬ 
files from our new Subaru data (Section 3), and examine the 
consistency of complementary shear and magnification mea¬ 
surements. In Section 4.2, we apply the 2D inversion method 
developed in Section 2 and reconstruct the projected 2D mass 
distribution from joint shear+magnification measurements. 

4.1. Weak-lensing Profiles of A1689 

A1689 exhibits a small offset d 0 g ~ 5 kpc/) -1 (~ 2.3") 
between the BCG and X-ray centroids (Table 1), ensuring a 
well-defined center. The X-ray and SZE centroids agree to 
within 1" (Table 1). Here we will adopt the BCG position as 
the cluster center for a radial profile analysis. 

r [/7 1 kpc] 


200 500 1000 2000 



Figure 5. Azimuthally averaged cluster weak-lensing profiles obtained from 
Subaru multi-color observations of A1689. The upper panel shows the tan¬ 
gential reduced shear profile (black squares) based on the full background 
sample. The lower panel shows the magnification-bias profile (red circles) 
of a z' - band limited sample of red background galaxies. For each observed 
profile, the shaded area represents the joint reconstruction ( 68 % CL) from 
the combined shear+magnification measurements. The horizontal bar (cyan 
shaded region) shows the constraints on the unlensed count normalization 
estimated from the source counts in cluster outskirts. 

We derive azimuthally averaged radial profiles of tangen¬ 
tial reduced shear ( g + ) and magnification bias {nf) from 
Subaru data. We calculate the lensing profiles in A| )itl = 
13 discrete radial bins, spanning the range [0 m in,$max] = 
[l',18'] with a constant logarithmic spacing, A Inf? = 

ln(0 max/^min)/-^bin — 0.22. The ilHlCr Tcldicll limit r m in — 

D\0mm ~ 129kpc/i -1 is sufficiently greater than the Ein¬ 
stein radius #Ein = 47.0" =b 1.2" (z s = 2; Table 1), and it 
also satisfies r m i n > 2d 0 ff — 10kpc h~ l , so that the miscen- 
tering effects on mass profile reconstructions are negligible 


r [h ^pc] 



Figure 6. Surface mass density profile E(0) (upper panel, red squares) de¬ 
rived from a Subaru ID weak-lensing analysis of the combination of shear 
and magnification measurements shown in Figure 5. The lower panel shows 
the corresponding cumulative mass profile M 2 d(< 0) (red squares). The 
gray area in each panel represents the best-fit projected Navarro-Frenk- 
White profile (68% CL) for the mass profile solution E(0). 

(Johnston et al. 2007; Umetsu et al. 2011a; Du & Fan 2014). 
The outer boundary 0 max = 18', or r max = A6» max ^ 
2.3 Mpc /i -1 , is large enough to encompass the entire virial 
region with r v ; r ~ 2 Mpc /A 1 (Umetsu & Broadhurst 2008), 
but sufficiently small compared to the size of the Suprime- 
Cam field of view so as to ensure accurate PSF anisotropy 
correction. The number of bins N\,i n = 13 is chosen such 
that the detection signal-to-noise ratio (S/N) is of the order of 
unity per bin, which is optimal for an inversion problem. 

In this work, we follow the prescription outlined in Section 
3.2.2 of Umetsu et al. (2014) to perform magnification mea¬ 
surements using the Subaru B /%:-^-selected red galaxy sam¬ 
ple (Table 4), which exhibits a clear depletion signal (Figure 
4). We have properly accounted and corrected for masking of 
background galaxies due to cluster galaxies, foreground ob¬ 
jects, and saturated pixels (see also Section 2.3.2). Unlike 
the nonlocal distortion signal, the magnification signal falls 
off sharply with increasing cluster radius. We thus estimate 
the count normalization and slope (n jL , a) from the source 
counts in cluster outskirts (Umetsu et al. 2011b, 2012, 2014; 
Medezinski et al. 2013), specifically at 12'(~ r2ooc) < 0 < 

f^max- 

Figure 5 shows the radial profiles of (g + ,n fl ). A clear de¬ 
pletion of red galaxies is seen toward the center owing to ge¬ 
ometric magnification of the sky area. The statistical signifi¬ 
cance of the detection of the tangential distortion is 22er. The 
detection significance of the magnification signal is 9er, which 
is ~ 40% of that of distortion. 

Here we construct the radial mass profile of A1689 from a 
joint likelihood analysis of shear and magnification measure¬ 
ments (Figure 5), using the method of Umetsu et al. (201 lb). 
We have 26 constraints {c/ +)i , n„ j}^ 11 in 13 log-spaced clus¬ 
tercentric radial bins. The model is described by N^in + 1 = 
14 parameters, {E min , , where E min = E(< d min ) 

is the average surface mass density interior to f? ttlin , and E, 
is the surface mass density averaged in the ith radial bin. To 
perform a reconstruction, we express the lensing observables 
(g + , p -1 ) in terms of E using the relations given in Appendix 
B. Additionally, we account for the calibration uncertainty in 
the observational parameters c = (( W) g , fw.g- (FF)^, n^,, a) 
as given in Tables 3 and 4. Following Umetsu et al. (2014), 
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we fix fw, g to the observed value (Table 3). 

The results are shown in Figures 5 and 6. The ML solu¬ 
tion has a reduced y 2 of 11.5 for 12 degrees of freedom (dof), 
indicating good consistency between the shear and magnifi¬ 
cation measurements having different potential systematics. 
This is demonstrated in Figure 5, which compares the ob¬ 
served lensing profiles with the respective joint reconstruc¬ 
tions (68% CL). The resulting mass profile 11(9) is shown in 
the upper panel of Figure 6. The error bars represent the ler 
errors from the diagonal part of the total covariance matrix 
C (Umetsu et al. 2014). The corresponding cumulative mass 
profile is shown in the lower panel of Figure 6. 

4.2. Weak-lensing Mapmaking ofA1689 

We apply our 2D inversion method (Section 2) to our new 
Subaru observations (Sections 3) for obtaining an unbiased 
recovery of the projected matter distribution E(0) in A1689. 
In this approach, we combine the observed spatial shear pat¬ 
tern (<7i(0), <72(0)) with the azimuthally averaged magnifica¬ 
tion measurements {n M !: }'^ b j" (Section 4.1), which impose 
a set of azimuthally integrated constraints on the underlying 
11(6) field, thus effectively breaking the mass-sheet degener¬ 
acy. The algorithm takes into account the nonlinear subcritical 
regime of the lensing properties. 

For mapmaking, we pixelize the lensing fields into a 56 x 56 
grid with Ad = 0.5' spacing, covering the central 28' x 28' 
field. The model m = (s , c) is specified by lV p i x = 56 2 pa¬ 
rameters, s = {£(0 m )}^L*, and a set of calibration param¬ 
eters c to marginalize over. We utilize the FFTW implemen¬ 
tation of fast Fourier transforms (FFTs) to calculate 7^ (6) 
from Koo(0) using Equation (6). To minimize spurious alias¬ 
ing effects from the periodic boundary condition, the maps 
are zero padded to twice the original length in each spatial 
dimension (e.g., Seljak 1998; Umetsu & Broadhurst 2008). 

We use a top-hat window of 9f = 0.4' (Section 2.3.1) to 
average over a local ensemble of galaxy ellipticities (N = 
Trrigdf ~ 10; Merten et al. 2014) at each grid point, ac¬ 
counting for the intrinsic ellipticity distribution of background 
sources. The filter size corresponds to an effective reso¬ 
lution of 2D\9{ ~ 100kpcfi -1 at the cluster redshift. To 
avoid potential systematic errors, we exclude from our analy¬ 
sis (Section 2.4.1) those pixels lying within central 9 cut = V 
where 11(9) can be close to or greater than the critical value 
£ c , as well as those containing no background galaxies with 
usable shape measurements. For distortion measurements 
(<7i(0), <72(0)) from the full background sample (Table 3), 
this leaves us with a total of 3093 usable measurement pix¬ 
els (blue points in Figure 7), corresponding to 6186 con¬ 
straints. For magnification measurements, we have 13 az¬ 
imuthally averaged constraints i j in log-spaced clus¬ 
tercentric annuli (Figure 7). The total number of constraints 
is thus IVdata = 6199, yielding A data - A pix = 3063 dof. 

In Figure 8, we show the resulting 11(9) field reconstructed 
from a joint analysis of the 2D shear and azimuthally av¬ 
eraged magnification data. The 7 2 value for the ML solu¬ 
tion is x 2 (rh) = 4046 for 3063dof. Here, for visualiza¬ 
tion purposes, the 11(9) field is smoothed with a Gaussian 
of FWHM = V. The main mass peak coincides well with 
the cluster center. The projected mass distribution is elon¬ 
gated in the north-south direction (Figure 1 ; see also Section 
7.1) and very similar to the distribution of cluster members 
(Kawaharada et al. 2010). 

In Figure 9, we compare the projected mass profiles 



x (arcmin) 


Figure 7. Spatial distribution of weak-lensing constraints averaged onto a 
grid of 56 x 56 pixels, covering a field of 28' x 28' centered on the BCG. Each 
point represents a single pixel with two-component reduced shear constraints 
( g \, ( 72 ) averaged within a top-hat region with radius Of = 0.4/. We exclude 
from our analysis those pixels lying within the inner 0 cu t = V region (red 
circle) and those having no background galaxies with usable shape measure¬ 
ments (see Figure 1). There are 3093 pixels with reduced-shear constraints, 
yielding 6186 constraints from 2D shear measurements. Azimuthally av¬ 
eraged magnification constraints are obtained in 13 logarithmically-spaced, 
clustercentric annuli spanning the range [0mm, 0max] = [1 / , 18 ; ]. 

11(9) obtained from our ID and 2D analyses of the 
shear+magnification data. Here we have used the method 
described in Appendix C to construct an optimally weighted 
radial projection of the £ map. Our ID- and 2D-based £ 
profiles are consistent within ler at all cluster radii, and both 
are in good agreement with the ID results of Umetsu et al. 
(2011b) from the joint shear+magnification analysis of the 
Subaru Vi' data. Similarly, our ID and 2D weak-lensing re¬ 
sults are in excellent agreement with each other in terms of 
the cumulative mass M- 2 d(< 9) as shown in Figure 10. 

5. HST STRONG-LENSING ANALYSIS 
5.1. Image Systems 

A1689 has been a subject of detailed strong-lensing studies 
by numerous authors (e.g., Broadhurst et al. 2005b; Halkola 
et al. 2006; Limousin et al. 2007; Coe et al. 2010; Diego et al. 
2015). Thus far, a total of 61 multiple-image candidate sys¬ 
tems of 165 images were identified from extremely deep opti¬ 
cal and near-infrared data from HST and Subaru (Diego et al. 
2015). 

To study global structural properties of the cluster, we fo¬ 
cus our strong-lensing analysis on the principal modes of the 
cluster mass distribution, responsible for the massive, smooth 
halo component (see Section 7.1.2). For this aim, we con¬ 
servatively select a subset of systems based on the following 
criteria: i) We use only spectroscopically confirmed systems. 
ii ) We consider only systems whose members were consis¬ 
tently identified in different studies, iii) We limit our analysis 
to those lying within 80"from the BCG, so that multiple im¬ 
ages spread fairly evenly over the analysis region, iv) We dis¬ 
card systems of very close pairs. They are primarily sensitive 
to substructures rather than the principal modes of the mass 
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Figure 8. Projected mass distribution E(0) of A1689 reconstructed from a 
Subaru weak-lensing analysis of 2D gravitational shear and azimuthally aver¬ 
aged magnification data. The mass maps is 2W x 28' in size (3.6 Mpc h~ 1 on 
a side) and centered on the BCG. The color bar indicates the lensing conver¬ 
gence k = (E(T 1 )E, scaled to the mean depth of weak-lensing observations, 
1/<E - ) = 4.66x 10 15 /iMq Mpc -2 . For visualization purposes, the mass 
map is smoothed with a 1' FWHM Gaussian. North is to the top, east to the 
left. The horizontal bar represents 1 Mpc h~ 1 at the cluster redshift. 

distribution, which we are interested in. 

These criteria leave us with 12 systems (ID 1, 2, 4, 5, 6, 7, 
11, 15, 18, 22, 24, 29, according to the original notation in 
Broadhurst et al. (2005b)), for a total of 44 multiple images 
spanning the range 1.4"-72.3" in cluster radius. 

5.2. PlXELENS Free-form Mass Reconstruction 

Free-form models describe the lens on a grid of pixels or a 
set of basis functions, allowing for a wide range of solutions 
(Coles 2008). We have performed a free-form strong-lensing 
analysis of the central region using the PixeLens software 
(Saha & Williams 2004), which produces pixelated maps of 
the surface mass density. Each map is constrained to exactly 
reproduce the positions and parities of all given multiple im¬ 
ages. PixeLens generates a statistical ensemble of models 
through which uncertainties and degeneracies in solutions can 
be explored (Coles 2008). 

Our PixeLens analysis procedure largely follows Sereno 
& Zitrin (2012) and Sereno et al. (2013). To determine robust 
sampling strategies optimized to recover the smooth cluster 
signal, we tested the PixeLens algorithm using simulated 
sets of multiple images in analytic lenses. The results suggest 
that the best strategy is to limit each analysis to three image 
systems, for a total of a dozen of images, and to reconstruct 
maps with ~ 10 pixels in the radial direction, avoiding over- 
sampling (Lubini & Coles 2012). We thus divide the strongly - 
lensed images in four groups of three systems each and ana¬ 
lyze each group separately. We end up with four triples con¬ 
sisting of systems 1, 5, and 11 (11 images), systems 2, 6, and 
22 (11 images), systems 4, 15, and 29 (12 images), and sys- 
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Figure 9. Comparison of projected mass density profiles E(r) derived 
from our Subaru ID weak-lensing analysis (squares; Section 4.1), Subaru 
2D weak-lensing analysis (circles; Section 4.2), and free-form strong-lensing 
analysis of HST data (triangles; Section 5). The cyan shaded area represents 
the mass profile with lcr uncertainty from a strong-lensing analysis of Broad¬ 
hurst et al. (2005b) based on the light-traces-mass (LTM) assumption. The 
gray shaded area shows the strong-lensing results (68% CL) from Diego et al. 
(2015) using a hybrid scheme combining both free-form grid and LTM sub¬ 
structure components. The diamonds with error bars show the results from 
our earlier ID weak-lensing analysis (Umetsu et al. 2011b) based on Sub¬ 
aru Vi' data. Good agreement between the strong and weak lensing results 
is seen in the region of overlap. There is also good agreement between the 
different lensing methods and data sets. 
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Figure 10. Comparison of projected cumulative mass profiles M 2 D (< r ) of 
A1689 derived from our Subaru ID weak-lensing analsyis (squares; Section 
4.1), Subaru 2D weak-lensing analysis (circles; Section 4.2), and HST strong- 
lensing analysis (triangles; Section 5). 


terns 7, 18, and 24 (10 images). Image systems with similar 
configurations are divided into different groups. 

For each group, we compute 500 k maps within 80" from 
the BCG on a circular grid of 349 pixels (10 pixels along 
the radial direction) with a pixel size of 8" (~ 17.2 kpc ft -1 ). 
These optimal settings allow us to avoid the known problem 
of too flat density profiles recovered with PixeLens model- 
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ing (see Grillo et al. 2010; Umetsu et al. 2012), which oth¬ 
erwise could bias cluster mass estimates. As discussed by 
Grillo et al. (2010, see their Appendix), this bias can arise 
from a combination of the mass-sheet degeneracy (Schneider 
& Seitz 1995) and the assumed prior on the positive definite¬ 
ness of every pixel of the surface mass density map. 

In the following, we restrict our analysis to the region where 
the cluster mass distribution is accurately recovered by PlXE- 
Lens. We exclude the central 20kpc h~ 1 region to minimize 
the effects of miscentering and baryonic physics (Umetsu 
et al. 2012, 2014). For each group of reconstruction, we 
determine the outer cutoff radius beyond which the logarith¬ 
mic density slope is steeper than -2, the asymptotic minimum 
slope for the projected Navarro-Frenk-White density profile 
(NFW, Navarro et al. 1997). The maximum radius is 63.7" 
(188 mass pixels) in three cases and 54.9" (140 mass pixels) 
for the group with the triple 4-15-29. 

5.3. Comparison of Weak and Strong Leasing Results 

We show in Figure 9 the radial mass distribution of A1689 
from our HST strong-lensing analysis. The results are shown 
along with the previous strong-lensing results by Broadhurst 
et al. (2005b) and Diego et al. (2015), as well as with indepen¬ 
dent weak-lensing results from shear and magnification infor¬ 
mation (Sections 4.1 and 4.2). The strong-lensing model of 
Broadhurst et al. (2005b) is based on the light-traces-mass 
(LTM) assumption, so that the HST photometry of cluster 
red-sequence galaxies was used as an initial guess for their 
lens solution. Diego et al. (2015) used a hybrid (free-form 
+ LTM) approach combining Gaussian pixel grid and clus¬ 
ter member components for describing large- and small-scale 
contributions to the deflection field, respectively. They con¬ 
strained the range of solutions with sufficient accuracy to al¬ 
low the detection of new counter images for further improving 
the lensing solution of A1689. This comparison shows clear 
consistency among a wide variety of lensing methods with 
different assumptions and potential systematics, demonstrat¬ 
ing the robustness of our results (see also Figure 10). Excel¬ 
lent agreement is also found between our strong-lensing mass 
profile and that of Limousin et al. (2007). 

6 . TRIAXIAL MODELING OF THE CLUSTER MATTER 
DISTRIBUTION 


(Navarro et al. 1996, 1997) 


p{r) 


Ps _ 

(r/r s )(l +r/r s ) 2 


(27) 


with p s the characteristic density and r s the inner characteris¬ 
tic radius at which the logarithmic slope of the density profile 
is -2. We generalize the spherical NFW model to obtain a tri- 
axial density profile by replacing r and r s with the respective 
ellipsoidal radii R and R s , defined such that 



Y 2 v 2 

^ + Z 2 , (28) 

<& q b 


where q a = a/c and q b = b/c (a < b < c) are the minor- 
major and intermediate-major axis ratios, respectively. ls The 
corresponding eccentricities are e a = — qf and e b = 

Vi — q 2 . The degree of triaxiality is defined as T = e 2 /e 2 
(Sereno et al. 2013). 

We define an ellipsoidal overdensity radius f? A (e.g., Cor- 
less et al. 2009; Sereno & Umetsu 2011; Buote & Humphrey 
2012b) such that the mean interior density contained within 
a ellipsoidal volume of semimajor axis f? A is A x p c . The 
total mass enclosed within f? A is Ma = (47 t/3) \q a q b p c R\. 
We use A = 200 to define the halo mass, M 2 ooc- The 
triaxial concentration parameter is defined by c 2 ooc = 
7?200 c/-Rs- The characteristic density is then expressed as 
Ps = M A /(^Ttq a q b R 3 A ) x <4/[ln(l + c A ) - c A /(l + c A )] 
(Buote & Humphrey 2012b). 

A triaxial halo is projected on to the sky plane as elliptical 
isodensity contours (Stark 1977), which can be expressed as 
a function of the intrinsic halo axis ratios (a/c, b/c) and ori¬ 
entation angles (i9, />, f ) with respect to the observer’s line of 
sight. Here we adopt the z-x-z convention of Euler angles to 
be consistent with Stark (1977) (see, e.g., Sereno et al. 2012). 
The angle 0 describes the inclination of the major (Z) axis 
with respect to the line of sight. 

For a given projection, the elliptical projected mass distri¬ 
bution can be described as a function of the elliptical radius / 
defined in terms of the observer’s coordinates (X', Y') in the 
plane of the sky: 

-| vv/2 

c 2 = -7 {jX' 2 + 2 kX'Y' + IY' 2 ) = (29) 

/ q\x q±Y 


Since we can only observe clusters in projection, determin¬ 
ing the intrinsic 3D shape and orientation of an aspherical 
cluster is an intrinsically underconstrained problem (Sereno 
2007). In this section, we describe the modeling of the 3D 
cluster matter distribution as an ellipsoidal halo following 
Sereno et al. (2013). In this approach, we exploit the com¬ 
bination of X-ray and SZE observations to constrain the elon¬ 
gation of the ICM along the line of sight. We use minimal ge¬ 
ometric assumptions about the matter and gas distributions to 
couple the constraints from lensing and SZE/X-ray data. The 
parameter space is explored in a Bayesian inference frame¬ 
work. This multi-probe method allows us to improve con¬ 
straints on the intrinsic shape and orientation of the cluster 
mass distribution without assuming HSE. 

6.1. Matter Distribution 

We model the cluster mass distribution with a triaxial NFW 
density profile as motivated by cosmological A'-body simula¬ 
tions (Jing & Suto 2002; Kasun & Evrard 2005). The radial 
dependence of the spherical NFW density profile is given by 


where q_ LX and q ±Y (q±x > Q±y) are 


q±x = 


2 

Y = 


2 / 


3+1- V / (j-0 2 + 4fc 2 ’ 
_2/_ 

j + 1 + VU ~ 0 2 + 4fc 2 ' 


(30) 


with 


J c 2 

j = cos" 1 t? ( cos 2 ( t ) + sin 2 < 
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C C , 2a 
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k = sin <b cos 6 cos i) 
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sm 


1 + 777 cos 
b z 
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(31) 


/ = sin 2 d 


77 sin 0+77 cos (f> + cos tL 
a z b z 


IS The intrinsic axis ratios (q a ,qb) here correspond to (tlDM.ai^DM b) 
of Limousin et al. (2013) in their notation. 
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Here we have chosen the new coordinate system (X" : Y") 
such that the X" axis is aligned with the major axis of the pro¬ 
jected ellipse. The minor-major axis ratio q± = q±y/q±x of 
the elliptical density contours is given by 19 


q±(a/c, b/c,id,(p) 


j + i-VU-i) 2 + 4fc 2 

j + l+ ^{j-l) 2 +±k 2 


(32) 

The principal axes of the isodensities are rotated by an angle 
tp with respect to the projection on to the sky of the intrinsic 
major axis Z, where 2 ip = arctan[2/c/(j — l )] (Sereno 2007). 
As observable parameters to describe the projected mass dis¬ 
tribution, we use the ellipticity 


e = 1 -q± (33) 

and the position angle i p t of the projected major axis. 

The projected surface mass density X(/) as a function of 
the elliptical radius ( is related to the triaxial density profile 
p(R) by (Stark 1977) 

2 f°° p{R)RdR 2 R s r°° p(R s x)xdx 

c _ v7 J„ 

(34) 

where £ = q±x( = \JdX" 2 -(- Y" 2 /q 2 ^ is the observable el¬ 
liptical radius, and £ s = q±xRs is the observable scale length 
(semi-major axis) in the sky plane (Sereno 2007). The quan¬ 
tity Z|| = -ff s /v/7 represents the line-of-sight half length of 
the ellipsoid of radius R = R s (Sereno 2007). It is useful to 
introduce the dimensionless scale factor ey that quantifies the 
extent of the cluster along the line of sight (Sereno 2007), 

e\\(a/c,b/c,&,<(>) = j- = f —) f~ 3/A . (35) 

Ss \9a9b/ 


The larger ey, the larger the elongation along the line of sight. 
The quantity ey corresponds to the inverse of the elongation 
parameter ex of Sereno (2007): ey = 1/eA- 
For a self-similar model p(R ) = p s fsD{R/Rs), the pro¬ 
jected mass density profile is expressed as 


2-R s p s f°° f 3 n(x)xdx 
~7T Jc/R, y/x 2 - (£/&) 2 


S s /2d(?/6)! 

(36) 


where we have defined the scale surface mass density 

£ s = 2 p s R s /s/f = 2p s £ s ey = 2/geoPsyTjX^s (37) 


with /geo = ey /y/q± (Sereno et al. 2010). Since r Sj2 D = 
^/q±^s is the geometric-mean scale radius in projection, 
the geometrical factor / geo represents the degree of cor¬ 
rection due to the line-of-sight elongation of the cluster. 
The halo mass, M 2 oo c , can then be expressed as M 2 00c = 
(47r/3)200p c (c 20 oc?’s,2D) 3 /geo. In this work, we employ 
the radial dependence of the projected NFW profile f-m(x) 
as given by Wright & Brainerd (2000). For / geo = 
1, this reduces to a projected (circular or elliptical) mass 
model. An elliptical mass density model can be described by 
(M 2 ooc, C 2000 e , ipe) (Oguri et al. 2010; Umetsu et al. 2012). 


6.2. Intracluster Gas 


19 Note the projected axis ratio q i is equivalent to 1/ep of Sereno (2007). 


Both observations and theory indicate that the ICM density 
is nearly constant on a family of concentric, coaxial ellipsoids 
(Kawahara 2010; Buote & Humphrey 2012a, b). Although 
modeling both the gas and matter distributions as ellipsoids 
with constant axis ratios is not strictly valid for halos in HSE 
(Sereno et al. 2013), an ellipsoidal approximation for the ICM 
is suitable when systems with modest eccentricities are con¬ 
sidered (Lee & Suto 2003). 

Following Sereno et al. (2013), we make a few simplifying 
but non-informative working hypotheses to relate the matter 
and gas distributions. First, we assume that the matter and 
gas distributions in the cluster are ellipsoidal with constant 
but different axis ratios and co-aligned with each other. Sec¬ 
ond, the two distributions are assumed to have the same de¬ 
gree of triaxiality, that is, T(q a ,qb) = T ICM (<7q CM , 9b° M ) 
with T icm = (e^/e 1 ™) 2 = [1 - (^ CM ) 2 ]/[1 - (gi CM ) 2 ] 
and (?^ CM < qJ CM . If two ellipsoids have the same de¬ 
gree of triaxiality, then the misalignment angle between their 
major axes in the plane of the sky is zero (Romanowsky & 
Kochanek 1998), which is consistent with what has been ob¬ 
served in A1689 (Sereno & Umetsu 2011; Sereno et al. 2012). 
If T = T icm , we have the following relation for the ratio of 
eccentricities between ICM and matter (Sereno et al. 2013): 

ef M /e a = ef M /e b = e ICM /e. (38) 

The intracluster gas in HSE is rounder than the underlying 
matter distribution: e ICM /e ~ 0.7 (Lee & Suto 2003). 

With these assumptions, the number of independent axis 
ratios is reduced to three. Here we use q a , q b , and q ^'' xl as free 
parameters. Hence, the intermediate-major axis ratio <?J CM of 
the ICM is determined by T(q a , 9b) and 9* CM : 


ICM 

9b 



(39) 


Finally, as supported by both theory and observations, we as¬ 
sume that the gas distribution is rounder than the matter dis¬ 
tribution: q a < 9), cm . 

Under these hypotheses, the projected matter and gas dis¬ 
tributions of the cluster have different ellipticities (e ^ e ICM ) 
and elongations (ey 7^ e|| CM ) but share the same orientation 
of the projected major axis, ip e = U) CM . There are a to¬ 
tal of six parameters ( q a , q b , 9 ^ CM , </, ip) needed to describe 

the intrinsic shape and orientation of the cluster system, com¬ 
pared to four observable geometric constraints, (e, e ICM , tp e = 
^ CM ,ef M ). 


6.3. Bayesian 3D Inversion 

In our analysis, the cluster model p is defined by seven fun¬ 
damental parameters describing the total matter ellipsoid and 
one parameter determining the shape of the ICM halo: 

P = (M 20 0c, C200c, 9a, 9b, id, <t>, ip, 9a CM )- (40) 

Hence, the overall ellipsoidal model has eight free parame¬ 
ters. On the other hand, 2D lensing constraints reduce to four 
parameters (Sereno & Umetsu 2011), (k s , £ s , e, ip e )- A joint 
X-ray and SZE analysis of the ICM yields two additional con¬ 
straints (Sereno et al. 2013), namely the ellipticity e ICVI of 
the ICM in projection and the elongation e|i CM of the ICM 
along the line of sight. Accordingly, combined lensing and 
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X-ray/SZE data sets effectively provide six observationally 
accessible parameters, 

o = (« s ,e s ,e,^,e ICM ,ef| CM ). (41) 

That is, the problem is underconstrained. 

To make robust inference on the intrinsic properties of the 
cluster, we use a forward modeling approach with Bayesian 
inference for this underconstrained inversion problem (Sereno 
et al. 2013). The observational parameters o = o(p) can be 
uniquely specified by the intrinsic parameters p. The total 
likelihood function of combined lensing and X-ray/SZE ob¬ 
servations can be formally written as (Sereno et al. 2013) 

C[o(p)\ = £ G l x £icm (42) 

with £ C i the likelihood function of lensing observables and 
£icm that of X-ray/SZE observables. 

6.4. Priors 

For our base model, we use uninformative priors for the 
intrinsic parameters p. We adopt flat priors of q m j n < q a < 1 
and q a < qb < 1 for the intrinsic axis ratios of the matter 
distribution, where g m ; n is introduced to exclude models with 
extremely small axis ratios because such configurations would 
be dynamically unstable and not expected for cluster halos. 
The probability functions can then be expressed as P(q a ) = 
1/(1 -g min ) for g min <q a < 1 and P(q b \q a ) = 1/(1 -g a ) _1 
for q b > q a ■ In what follows, we fix q min = 0.1 (Oguri et al. 
2005; Sereno et al. 2013). Alternatively, we may consider 
the axis-ratio priors that follow distributions obtained from 
ACDM A' -body simulations (Jing & Suto 2002). 

For the minor-major axis ratio of the ICM, we use a uni¬ 
form distribution in the interval q a < g* CM < 1 (see Sec¬ 
tion 6.2). The prior of -P(ga CM |g a ) can then be de¬ 

fined in a similar way to that of q b . For the orientation an¬ 
gles, we consider a population of randomly oriented halos 
with P( cost?) = 1 for 0 < cost? < 1 and P{(j>) = 1 /7r 
for —7 t/2 < (!) < </>/2. Finally, we employ uniform priors for 
the remaining parameters. 

7. MULTI-PROBE ANALYSIS OF A1689 

Here we apply the Bayesian inversion method outlined in 
Section 6 to our multiwavelength observations of A1689. The 
results are discussed in Section 8. 

7.1. Weak and Strong Lensing 

A full 2D lensing analysis is crucial for comparison with 
predictions of the properties of aspherical clusters (Oguri et al. 
2005). In this work, we have employed free-form methods for 
both weak- and strong-lensing mass reconstructions (Sections 
4 and 5), which provide a pixelated X map and its covariance 
matrix in each regime. 

In this subsection, we derive constraints on the projected 
halo properties (Section 6.1) from lensing data. We model 
the observed X field with a projected ellipsoidal NFW pro¬ 
file (Section 6.1), specified by (k S) £ S) e, tjt e ). Additionally, we 
include the halo centroid Q c as parameters to conservatively 
account for the degree of miscentering. 

7.1.1. Weak-lensing Data 



k s [z s =2.0] 

Figure 11. Marginalized posterior distribution for the projected NFW pa¬ 
rameters (k s , 0 s ) obtained from three different lensing data sets (see Table 
5) , namely weak-lensing-only (black; WL), strong-lensing-only (blue; SL), 
and combined weak and strong lensing (red shaded; WL+SL). For each case, 
the contour levels are at exp(—2.3/2) and exp(—11.8/2) of the maximum, 
corresponding to the lcr and 3cr confidence levels, respectively, for a Gaus¬ 
sian distribution. The scale convergence k s = E s /E c is normalized to a 
fiducial source redshift of z s = 2. 

The x' 2 function for the Subaru weak-lensing observations 
is expressed as (Oguri et al. 2005) 


N pix 

XwL = E n0m)-t(0 m ) {C~ X ) mn X(0 n )-E(0 n ) 


(43) 

where X = {X(0 m )} m X(j is the mass map from the 2D 
weak-lensing analysis (Section 4.2), C -1 is the inverse 
of the error covariance matrix, and the hat symbol de¬ 
notes a modeled quantity. The corresponding likelihood is 

£wL(K s ,( sl e,i,0c) oc exp(-Xw L / 2 )- 

Figure 11 shows the results in terms of the marginalized 
posterior distribution for the scale convergence, n s = X s /X c , 
and the scale radius, 0 S = £ S /D\. Table 5 summarizes 
marginalized constraints on the individual parameters. In the 
present study, we employ the robust biweight estimators of 
Beers et al. (1990) for the central location (mean) and scale 
(standard deviation) of the marginalized posterior distribu¬ 
tions (e.g., Sereno & Umetsu 2011; Umetsu et al. 2014). 


7.1.2. Strong-lensing Data 

Mass maps derived from strong lensing exhibit a high de¬ 
gree of correlation between adjacent regions. The problem 
is exacerbated for parametric methods, which model the total 
mass distribution by a superposition of lens components as¬ 
suming parametric density profiles. This also persists in free¬ 
form modeling (Lubini et al. 2014), albeit to a lesser degree. 

The degree of correlation can be examined by an eigenvalue 
analysis. Let us decompose the C matrix as C = UAU _1 , 
with A the diagonal matrix of eigenvalues and U the unitary 
matrix of eigenvectors. The first few eigenvalues describe the 
principal modes of variation of the mass model (Lubini et al. 
2014; Mohammed et al. 2014). Large eigenvalues correspond 
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Table 5 

Parameters of the projected NFW model constrained from lensing observations 


Data a 

K s b 

(') 

e d 

V-e e 

(deg) 

0c f 

<") 

WL 

0.97 ± 0.16 

1.74 ± 0.27 

0.29 ±0.07 

14.2 ±8.4 

-1.2 ±3.0, 4.9 ±4.1 

SL 

0.73 ±0.14 

3.00 ±0.90 

0.27 ± 0.09 

13.0 ±9.8 

-0.8 ±1.9,-4.8 ±2.3 

GL 

1.03 ±0.11 

1.70 ±0.20 

0.29 ±0.05 

11.4 ± 4.9 

0.0 ± 1.3, —1.9 ± 1.4 


a WL: weak lensing shear and magnification; SL: strong lensing; GL: combined strong lensing, weak- 
lensing shear and magnification. 

b Scale convergence, k, s = E s /E c , normalized to a reference source redshift of z s = 2. 
c Projected scale radius of the elliptical NFW model measured along the major axis. 
d Projected mass ellipticity, e = 1 — q±, with q± the projected minor-major axis ratio. 
e Position angle of the major axis measured east of north. 
f Halo centroid position relative to the BCG position. 


to massive pixels, namely, those composing the inner part of 
the mass distribution that is best constrained by strong lens¬ 
ing. The ordered list of eigenvalues progressively decreases 
with increasing rank and drops abruptly near the maximum 
rank, indicating a high degree of correlation (Figure 12). 



Figure 12. Ordered eigenvalues A of the covariance matrix for the PlXE- 
LENS mass reconstruction. The vertical red line indicates the maximum rank 
considered for our analysis, AT ma x = 2Ni m , i.e., the number of observa¬ 
tional constraints on the image position. The blue horizontal line shows the 
minimum k 2 value found in the ensemble-averaged pixelated model. The 
green horizontal line shows the minimum k 2 value from the entire statistical 
ensemble of models generated by PixeLens. The results are shown for the 
covariance matrix as constrained by the systems 1, 5, and 11. 


Here, we employ a regularization approach to conserva¬ 
tively account for the high degree of correlation of the covari¬ 
ance matrix. This was first proposed by Umetsu et al. (2012) 
for the ID analysis of strong-lensing mass profiles. If the co- 
variance matrix C is not degenerate, we can construct a x 2 
function for each group of multiple images as 


XSL,Q! — [ 

m,n 

= £ 


V — V 


(c- 1 ) 

V / 7 




[£t/)m - (^u)r. 


(44) 




where £ m is the observed £ value of the mth pixel, (= 
Y2i Umi'Ei is the projection onto the eigenbasis, a runs over 
the four groups of images (Section 5), and the hat symbol is 
used to denote a modeled quantity. Each group has its own £, 
C, U, and A. Here we drop the index a on the right hand side 
to simplify the notation. 


In this approach, we limit ourselves to the principal modes 
and truncate the summation at iV max largest eigenvalues as 


N m 


X'SL.c 


£ 


[Pu)ro “ (^u)r 
A m 


(45) 


A natural choice for 7V max is the number of observational con¬ 
straints. We thus set N max = 2 N- lm with N- lln the number of 
multiple images used. The total \ 2 is given by 


XSL — £xSL,a- (46) 

a. 


We find that the eigenvalues before the drop range ap¬ 
proximately between the minimum re 2 value in the ensemble- 
averaged pixelated model and that found from the whole en¬ 
semble of models generated by PixeLens (Section 5.2). This 
is demonstrated in Figure 12. The 2iVj m -th eigenvalue lies ap¬ 
proximately in the middle of this range and sets a conservative 
scale. We checked the reliability and performance of this reg¬ 
ularization method using analytical models. 

Some multiple image systems share very similar configura¬ 
tions (e.g., systems 1 and 2). Such a redundancy is valuable 
for determining cosmological parameters (Lubini et al. 2014), 
or for improving the sensitivity to local substructures. As¬ 
signing a full weight to systems having similar configurations 
would inflate the relative contribution of strong lensing with 
respect to weak lensing. To avoid this, we multiply %g L by a 
weighting factor wsl, defined as the inverse of the geometri¬ 
cal average of the number of such redundant image systems. 
We find wsl = 2/3 for our analysis. The likelihood is then 
defined as £ S l(k s , £ s , e, ipt, 0 C ) ex exp(-u;sLXsL/ 2 )- 

The results are summarized in Table 5 and Figure 11. 


7.1.3. Combining Weak and Strong Lensing 

We now combine the weak- and strong-lensing likelihoods 
constructed in Sections 7.1.1 and 7.1.2, respectively, to jointly 
constrain the projected NFW parameters. The likelihood 
function £gl for the combined weak plus strong lensing data 
can be written as (Sereno & Umetsu 201 1) 

£ql = £wl x £ sl ex exp[— (xwl + w slXsl)/ 2 ]> (47) 

where Xwl an d Xsl are defined by Equations (43) and (46), 
respectively. 

Figure 11 shows that the scale radius ( 9 S ) and the scale 
convergence (n s ) are highly degenerate and anti-correlated. 
In particular, the scale radius is poorly constrained by strong 
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Table 6 

Integrated Comptonization Y 
parameter measured interior to a 
cylinder of radius r. 


Instrument 

r 

0 

Y{< r) 
(10 —10 sr) 

BIMA/OVRO 

1.5 

1.00 ±0.28 

BIMA/OVRO 

3.0 

2.64 ± 0.97 

SZA 

1.5 

1.11 ± 0.10 

SZA 

3.0 

2.83 ±0.42 

SZA 

4.5 

4.33 ±0.81 

SZA 

6.0 

5.50 ± 1.18 


lensing alone because of the limited coverage of multiple im¬ 
ages, 9 4S 1.1' (Section 5.2). The allowed range of 9 S lies 
well outside the region where the multiple images are ob¬ 
served. Thus, the inference of parameters by strong lens¬ 
ing requires an extrapolation well beyond the observed re¬ 
gion. For this reason, in the present study, we do not con¬ 
sider strong-lensing-only triaxial modeling (see Table 7). On 
the other hand, since the posterior distributions from the inde¬ 
pendent weak-lensing and strong-lensing analyses are com¬ 
patible, combining weak lensing with strong lensing provides 
improved parameter constraints (Table 5). 


7.2. Combined X-ray plus SZE Analysis 

With a known halo geometry (e.g., sphericity) and under 
the ideal gas assumption, the thermodynamic quantities of the 
ICM are overconstrained by X-ray and SZE data. This is be¬ 
cause the thermal pressure can be independently determined 
from thermal SZE data and X-ray spectroscopy/imaging data. 
We can therefore relax the assumption of spherical sym¬ 
metry to solve for the elongation of the ICM distribution 
(Sereno et al. 2012). Combining gravitational lensing and 
X-ray/SZE observations with minimal geometric assumptions 
(Section 6.2) allows us to break the degeneracy between mass 
and elongation for the total matter distribution (Sereno et al. 
2013). Such a multi-probe approach based on lensing and X- 
ray/SZE data is free from the assumption of HSE, compared 
to the lensing plus X-ray analysis, which relies on equilibrium 
assumptions between the gravitational potential and pressure 
components (see Limousin et al. 2013). 

In our multi-probe approach, the ICM distribution is mod¬ 
elled with an ellipsoidal parametric profile which can fit X- 
ray surface-brightness (Sx) and temperature (T x ) distribu¬ 
tions. Comparison with the SZE amplitude then determines 
the elongation ej| CM For example, for an isothermal plasma 
(De Filippis et al. 2005), we have 


l/ef M 


oc D\ 


Sx rj 

AlfzE A.y 


(48) 


with ATsze the SZE temperature decrement and A y the X- 
ray cooling function of the ICM. In this work, we rely on 
the X-ray data to constrain the ICM morphology in projection 
space; we use aperture-integrated constraints on the SZE sig¬ 
nal (Table 6) to determine the line-of-sight elongation ejj CM . 

Our X-ray data are taken from Sereno et al. (2012), who 
performed an X-ray analysis on Chandra and XMM-Newton 
observations. Here we briefly summarize essential results 
needed for this study. For details, we refer to Sereno et al. 
(2012). Sereno et al. (2012) showed that exposure corrected 
and point-source removed Chandra X-ray images in the 0.7- 
2.0 keV band are well described by concentric ellipses with 


ellipticity e x = 0.15 ± 0.03 and orientation angle ij) X = 
(12 ± 3) degrees measured east of north. Following Sereno 
et al. (2012) and Sereno et al. (2013), we model the 3D elec¬ 
tron density in the intrinsic coordinate system with the follow¬ 
ing parametric form (Vikhlinin et al. 2006; Ettori et al. 2009): 


n e 


n 0 



(49) 


where no is the central electron density, r c is the ellipsoidal 
core radius, r t (> r c ) is the ellipsoidal truncation radius, fi 
is the slope in the intermediate density regions, and 7 is the 
outer slope. The 3D gas density is parametrized as (Sereno 
etal. 2013) 


T = 


T 0 

1 + (R/tt) 2 


(50) 


where Tq is the central gas temperature, and t j t describes a 
temperature decline at large cluster radii. The parametriza- 
tions of Equations (49) and (50) were motivated by the ab¬ 
sence of cool-core features in our data. For further justifica¬ 
tion, see Section 5 of Sereno et al. (2012). 



Figure 13. Marginalized posterior probability distribution of the elongation 
e( CM as derived from the combined X-ray plus SZE analysis (Section 7.2). 


The thermal SZE provides a complementary measure of 
the thermal energy content in a cluster. In this study, 
we perform a self-consistent multi-scale analysis of high- 
significance 30 GHz interferometric SZE observations of 
A1689 obtained with the Berkeley-Illinois-Maryland Array 
(BIMA), the Owens Valley Radio Observatory (OVRO), and 
the Sunyaev-Zel’dovich Array (SZA). The BIMA and OVRO 
observations of A1689 are presented in LaRoque et al. (2006), 
while the SZA observations of A1689 are presented in Gralla 
et al. (2011). Owing to the different scales probed by the in¬ 
struments, we fit the OVRO/BIMA and SZA data separately 
using the spherical Arnaud et al. (2010) pressure profile. This 
profile is an adaptation of the generalized NFW pressure pro¬ 
file first proposed by Nagai et al. (2007), and first fitted to 
SZE observations in Mroczkowski et al. (2009). A joint fit 
to the OVRO, BIMA, and SZA data was also performed to 
determine the best-fit SZE centroid reported in Table 1. 

As in Mroczkowski et al. (2009), a model for the cluster and 
contaminating radio sources is computed in the image plane. 
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then Fourier transformed for comparison to the interferomet¬ 
ric data. The best-fit model and lcr confidence intervals are 
determined using a Markov chain Monte Carlo (MCMC) pro¬ 
cedure. The OVRO and BIMA data measure radial scales 
from 0.5'-4', while the SZA data probe radial scales from 
1'—6'. Bonamente et al. (2012) showed that the adoption of 
the Arnaud et al. (2010) profile versus other non-isothermal 
pressure profiles accurate out to r §ooc does not significantly 
impact the parameters derived from the fits when the radii for 
which the results are computed are at scales accessible to the 
instruments. 

A summary of the SZE data used is given in Table 6. The 
integrated Comptonization parameter Y(< r ) interior to a 
cylinder of radius r is written in terms of the electron den¬ 
sity and temperature profiles (Equations (49) and (50)) as 

Y= arkB f dQ f dlrieT (5i) 

m e c z J Qr J 

with ctt the Thomson cross section, ky>, the Boltzmann con¬ 
stant, m e the electron mass, and c the speed of light in vac¬ 
uum; Q r is the solid angle of the integration aperture. 

The model profiles given by Equations (49), (50), and 
(51) are then compared with combined X-ray surface bright¬ 
ness (Sx), X-ray spectroscopic temperature (Tx), and ther¬ 
mal SZE decrement (Y) observations. Briefly summarizing, 
the X-ray surface brightness profile observed by 

Chandra was extracted from Ns = 68 elliptical annuli out 
to an elliptical radius of £ = 900 kpc hy 0 (~ 5'), and the 
XMM-Newton temperature profile {Tx <t }fYy was measured 
in Nt = 5 elliptical annual bins out to £ = 900 kpc h? 0 
(Sereno et al. 2012). Thanks to the improved SZE analy¬ 
sis, the Y parameter is measured at several apertures from 
BIMA/OVRO and SZA data as summarized in Table 6. We 
find good consistency between the BIMA/OVRO and SZA re¬ 
sults at r = 1.5' and 3' where these independent data overlap. 
At an integration radius of r = 3', our results are also in ex¬ 
cellent agreement with Y(< 3') = (2.5 ± 0.6) x 10~ ln sr 
from 94 GHz interferometric observations with the 7-element 
AMiBA (Umetsu et al. 2009, their Table 5). 

The X-ray part of the \ 2 function can be written as (Sereno 
et al. 2012) 


N s 


Xx 


= £ 


Sx,i — Sx,i 

a S,i 


Nt / rp rft 

_|_ | -*-X,i -LX,i 

h v ^ 


2 

(52) 


with ( Sx , 7’y) model predictions for the corresponding X-ray 
observables and (as, or) their corresponding errors. 

The x 2 function for the SZE observations is written as 



A Yjj - AY A 
a X,ji J 


(53) 


where AY rt is the differential Y parameter for the yth instru¬ 
ment (BIMA/OVRO or SZA) in the /tli annular ring, AYj, = 
Yj(< ri+ 1) — Yj (< n), and ctaj-i is its lcr uncertainty. The 
Y values are sampled at every 1.5' (Table 6), which is suffi¬ 
ciently larger than the synthesized beam. Hence, differential 
AYmeasurements in adjacent annuli are approximately un¬ 
correlated given the annulus size considered. 

A combined analysis of the X-ray and SZE data is per¬ 
formed using the combined function x ‘ 2 = Xx + x|ze- The 


parameter space is explored using an MCMC approach as de¬ 
scribed in Sereno et al. (2013). Since parameter constraints 
on the n e and T models are dominated by the Chandra sur¬ 
face brightness and XMM-Newton temperature data, respec¬ 
tively, we find our results are fully consistent with those of 
Sereno et al. (2012) based on the same X-ray data. The best- 
fit central temperature (To = 9.8 ± 0.2keV, Sereno et al. 
2012) is in good agreement with the Suzaku X-ray results of 
Kawaharada et al. (2010). On the other hand, using the im¬ 
proved SZE data, we obtain tighter constraints on the elon¬ 
gation e|| CM . The resulting posterior distribution of e|| CM is 
shown in Figure 13. The posterior mean and standard devia¬ 
tion are ef,™ = 1.70 ±0.29. 


7.3. Multi-probe Deprojection 

Here we perform joint likelihood analyses of combined 
lensing and X-ray/SZE data, using different combinations of 
lensing data sets (Section 7.1). 

The likelihood £icm of the X-ray/SZE data is written 
in terms of two observable ICM parameters (Section 6.3), 
namely the ellipticity e ICM and line-of-sight elongation e|| CM 
of the ICM. Following Sereno et al. (2012, 2013), we include 
a nuisance parameter Aej) ys that quantifies the additional un¬ 
certainty on ej| CM , accounting for potential calibration sys- 
tematics in the X-ray/SZE measurements. It is assumed to 
follow a normal distribution with zero mean and standard 
deviation erjj ys = 0.07. Since the systematic uncertainty is 
quite small compared to the width of the marginalized pos¬ 
terior distribution P(ejj CM ) (Figure 13), the impact on the fi¬ 
nal results is minor. The X-ray/SZE part of the likelihood 
AcM(e ICM > e [| CM ; ^ e y yS ) is written as (Sereno et al. 2013) 


£icm = 


exp 


r e X _ e ICM)2 

^ a t,X 


\Fhto^x 
x P(ef, CM - Ae?, ys ) 


(54) 


■s/^oT 


exp 


i^TY 


sys 


where e x and a e x are the measured value of the ICM ellip¬ 
ticity and its uncertainty, respectively (Section 7.2). 

To perform a joint analysis with the X-ray/SZE data, we 
consider three different likelihood functions for the lensing 
part, namely, £ W l> £sl> and £gl = Ayl^sl, which are 
all functions of the projected NFW parameters K a ,^ s ,e,ijj e , 
and 9 C . Following Sereno et al. (2013), we exploit constraints 
from the X-ray analysis about the gas centroid 6 X and posi¬ 
tion angle ip x (Section 7.2), which are used as priors for the 
centroid 6 C and position angle ti; ( of the underlying halo (see 
Section 4 of Sereno et al. 2013). These priors are consistent 
with the geometric assumptions we have made in Section 6.2. 

For our base model, we use flat priors for the intrinsic axis 
ratios of the underlying halo (Section 6.4). We also consider 
an alternative prior distribution predicted by cosmological N- 
body simulations of ling & Suto (2002). For details, we refer 
to Sereno & Umetsu (2011) and Sereno et al. (2013). 

8. RESULTS AND DISCUSSIONS 

The resulting constraints on the intrinsic parameters for the 
underlying halo (M 20 oc, c 2 ooc, q a , Qb, cost?) are given in Ta- 
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Table 7 

Intrinsic parameters of the total matter distribution obtained using different data sets and different priors 


Data a 

Prior 

M 200 C 

(1O 15 A7 0 /i -1 ) 

C200c 

9a 

96 

COS$ b 

WL 

Spherical 

1.31 ± 0.11 

8.87 ± 1.11 

1 

1 

— 

WL 

Flat 

1.28 ±0.26 

10.70 ± 2.85 

0.39 ±0.18 

0.77 ± 0.15 

0.54 ± 0.29 

WL 

7V-body 

1.22 ±0.23 

9.15 ± 1.77 

0.47 ± 0.08 

0.66 ± 0.12 

0.60 ±0.30 

SL 

Spherical 

1.79 ±0.31 

8.69 ± 1.26 

1 

1 

— 

GL 

Spherical 

1.32 ±0.09 

10.10 ±0.82 

1 

1 

— 

GL 

Flat 

1.49 ±0.25 

10.30 ± 2.52 

0.45 ±0.20 

0.77 ± 0.14 

0.47 ±0.29 

GL 

W-body 

1.41 ± 0.19 

9.65 ± 1.54 

0.47 ± 0.08 

0.66 ± 0.12 

0.60 ±0.29 

WL + X/SZ 

Flat 

1.21 ± 0.19 

7.91 ± 1,41 

0.39 ±0.16 

0.56 ± 0.20 

0.93 ± 0.06 

WL + X/SZ 

W-body 

1.16 ±0.17 

7.42 ± 1.21 

0.40 ±0.08 

0.52 ± 0.12 

0.94 ± 0.05 

GL + X/SZ 

Flat 

1.24 ± 0.16 

8.36 ± 1.27 

0.39 ±0.15 

0.57 ± 0.19 

0.93 ± 0.06 

GL + X/SZ 

7V-body 

1.20 ±0.13 

7.89 ± 0.96 

0.40 ±0.08 

0.52 ± 0.12 

0.94 ± 0.05 


Note. — Intrinsic parameters of the total matter distribution of A1689 derived from a triaxial analysis of 
multiwavelength data sets, using spherical, flat, and JV-body priors on the distribution of axis ratios (g a , qt>). 


a WL: weak-lensing shear and magnification; SL: strong lensing; GL: combined strong lensing, weak-lensing 
shear and magnification; X/SZ: combined X-ray and SZE measurements. 
b Cosine of the angle between the major axis and the line of sight. 


Table 8 

Intrinsic shapes of the ICM distribution 


Priors 

—[CM 

"a 

„ICM 

% 

gtCM/e 

Flat 

0.60 ±0.14 

0.70 ±0.16 

0.87 ±0.07 

7V-body 

0.58 ±0.10 

0.65 ± 0.11 

0.89 ± 0.06 


Note. — Constraints on the intrinsic axis ratios 
l7a CM i9{, CM ) °f the ICM distribution and the relation 
with the total matter distribution (e ICM /e), obtained 
from the full triaxial analysis of combined weak/strong- 
lensing and X-ray/SZE data sets (Section 7). and 

e lCM/ e are derived parameters. 

ble 7, for different combinations of data sets and three dif¬ 
ferent priors on the axis-ratio distribution: (1) spherical prior 
(q a = Qb = 1); (2) flat distribution of axis ratios and random 
distribution of halo orientations (Section 6.4); (3) 7V-body 
ACDM predictions (Jing & Suto 2002). The baseline results 
for the combined weak/strong-lensing and X-ray/SZE analy¬ 
sis obtained with flat priors are shown in Figure 14. Table 8 
gives a summary of our baseline constraints on the intrinsic 
axis ratios of the ICM halo, (q l ^ M , q\ CM ), and on the ICM- 
to-matter ratio of halo eccentricities, e ICM /e. Table 9 lists the 
published (M 2 ooo C 200 c) measurements for A1689 based on 
the combination of both weak and strong lensing. For pre¬ 
vious compilations, see Comerford & Natarajan (2007, their 
Table Al) Fimousin et al. (2007, their Table 4), Umetsu & 
Broadhurst (2008, their Table 5), Corless et al. (2009, their 
Table 4), and Coe et al. (2010, their Table 2). 

8.1. Mass and Concentration 

8.1.1. Spherical Modeling 

The degree of concentration of A1689 has been a subject 
of controversy. Here we first compare the results obtained 
assuming a spherical NFW halo (Table 7) to those of pre¬ 
vious work. Our full 2D weak-lensing analysis based on 
Subaru BVRqi' z' data yields a projected concentration of 
C 200 c = 8.9 ± 1.1 (c v i r = 11.2 ± 1.4) at M 2 oo c = (1-31 ± 
0.11) x 10 15 Mq/i 1 . This is in excellent agreement with, 
and improved from, our earlier weak-lensing work: c 2 oo c = 
10.7^2 7 (Umetsu & Broadhurst 2008) and c 2 ooc = 10.2^2 q 
(U metsu et al. 201 lb), both of which are based on the joint 
analysis of shear and magnification data from Subaru Vi' 


imaging. 2,1 This accurate agreement comes in spite of using 
different data reduction procedures and mass reconstruction 
methods (Sections 2, 3, and 4). 

Combining weak and strong lensing reduces the uncer¬ 
tainties on the concentration. The HST strong-lensing data 
alone also favor a high degree of projected concentration, 
c 2 ooc = 8.69 ± 1.26, but with a somewhat higher halo mass, 
M 20 oc = (1-79 ± 0.31) x lO 15 M 0 /i _1 . The combined 
weak and strong lensing data yield c 2 oo c = 10.10 ± 0.82 at 
M 2 ooc = (1-32 ± 0.09) x 10 15 M 0 /i _1 , corresponding to the 
Einstein radius of (/j.-j,, = 52+®" at z s = 2. Our analysis thus 
reproduces the correct size of the observed Einstein radius 
(Table 1). These results are in good agreement with those of 
Umetsu & Broadhurst (2008) and Coe et al. (2010) (Table 9), 
in spite of using completely independent approaches to strong 
lens modeling (Section 5). Most recent weak-and-strong lens¬ 
ing studies of A1689 appear to converge toward c 2 oo c ~ 9-10 
with a typical measurement uncertainty of 10% (Table 9; with 
the spherical prior), thanks to the advanced analysis methods 
and greatly improved quality of data. 

8.1.2. Triaxial Modeling 

Including triaxiality weakens parameter constraints from 
lensing data (Oguri et al. 2005; Corless et al. 2009), compared 
to those derived assuming spherical symmetry. The parameter 
constraints become more degenerate and less restrictive be¬ 
cause of the lack of information of the halo elongation along 
the line of sight (Table 7). These trends are also found in the 
posterior distributions from our data (Tables 7 and 9). 

Now we consider the results from full triaxial analyses 
combining lensing with X-ray/SZE data. Table 7 shows that 
our posterior inference of the intrinsic parameters is insen¬ 
sitive to the assumed choice of priors (“Flat” or “A'-body”) 
when the line-of-sight information from X-ray/SZE data is 
combined with lensing, suggesting that the posterior con¬ 
straints are dominated by the likelihood (i.e., information 
from data) rather than the prior (Sereno et al. 2013). What¬ 
ever the assumptions regarding the axis ratios, we find the 

20 Umetsu & Broadhurst (2008) derived a n(6) map for the cluster us¬ 
ing an entropy-regularized maximum-likelihood combination of 2D shear 
and magnification maps. Umetsu et al. (2011b) derived a k(6) profile from 
a joint likelihood analysis of azimuthally-averaged shear and magnification 
measurements. 
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Figure 14. Marginalized posterior distributions for the intrinsic parameters of the triaxial cluster model obtained from a joint analysis of the weak/strong-lensing 
and X-ray/SZE data. In each panel, the contours levels are shown at exp(—2.3/2), exp(—6.17/2), and exp(—11.8/2) of the maximum, corresponding to the 
lcr, 2cr, and 3cr confidence levels, respectively, for a Gaussian distribution. In the qf, versus q a plane, the green solid (diagonal) and dashed (horizontal) lines 
represent prolate (q a = qb) and oblate (qj = 1) configurations, respectively, and the thick red line shows the lcr contour for the axis-ratio distribution in ACDM 
iV-body simulations of Jing & Suto (2002). 


posteriors (Table 7) to be statistically compatible with the 
predicted distribution c(M) for the full population of ha¬ 
los in ACDM cosmological simulations (Bhattacharya et al. 
2013; Meneghetti et al. 2014; Diemer & Kravtsov 2015). 21 
This is demonstrated in Figure 15 for the weak-lensing plus 
X-ray/SZE analysis and for the weak/strong-lensing plus X- 
ray/SZE analysis, both based on the uninformative priors. 
Here we adopt the median c-M relation obtained by Diemer 
& Kravtsov (2015) as a reference model for comparison. 

A1689 appears to be a high mass cluster of M 200 C ~ 
10 15 Mq/i _1 in the high-concentration tail of the predicted 

- 1 The theoretical predictions from Bhattacharya et al. (2013) and Diemer 
& Kravtsov (2015) are based on DM-only simulations, and those from 
Meneghetti et al. (2014) are based on nonradiative simulations of DM and 
baryons. 


c(M) distribution (Figure 15). The posterior tail at lower con¬ 
centrations of A1689 is only £ lcr away from the predicted 
median concentration (log 10 C 200 C — 0.58 ± 0.16; Figure 15). 
Our results are also in agreement with those obtained by a 
multi-probe analysis of Sereno et al. (2013) (see Table 9), who 
developed the triaxial inversion algorithm used in this work. 

The halo concentration and orientation are strongly corre¬ 
lated (Sereno & Umetsu 2011; Sereno et al. 2013). For the 
posterior range 0° < 1 / < 5° assuming a nearly perfect align¬ 
ment between the halo major axis and the line of sight, we 
find C 200 C = 7.4 ± 1.0 (6.7 ± 1.1) from weak/strong lensing 
(weak lensing) combined with the X-ray/SZE data. 


8.2. Intrinsic Shape and Orientation ofA1689 
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Table 9 

Published mass and concentration measurements of A1689 from combined weak and strong lensing 


Author 

M 200 C 

C200c 

Prior* 

External data b 


(to 15 tf 0 r 1 ) 




Spherical modeling 
Broadhurst et al. (2005a) 

1.20 ±0.13 

io-9ioi 

Spherical 


Halkola et al. (2006) 

1.58 ±0.14 

7.6 ± 0.5 

Spherical 

— 

Umetsu & Broadhurst (2008) c 

1.30 ±0.11 

10-l±o ® ± 2.2 

Spherical 

— 

Coe et al. (2010) 

This work 

| q+0.3 

0.2 

1.32 lb 0.09 

9.2 ± 1.2 
10.10 ±0.82 

Spherical 

Spherical 

— 

Triaxial modeling 

Oguri et al. (2005) d 

Sereno & Umetsu (2011) 

1 14 + 0- 26 

1.07 zb 0.23 

13 6 +1,8 
rj-O-10.5 

9.3 ±2.0 

Flat 

Flat 

— 

This work 

1.49 ±0.25 

10.30 ±2.52 

Flat 

— 

With line-of-sight information 





Corless et al. (2009) 

0.83 ±0.16 

12.2 ±6.7 

Flat + cos $ 

— 

Sereno & Umetsu (2011) 

0.99 ±0.17 

7.7 ± 1.1 

Flat + cos $ 

— 

Morandi et al. (2011) 
Sereno et al. (2013) 

This work 

1.81 ± 0.06 
0.93 ±0.12 
1.24 ± 0.16 

5.71 ± 0.47 

7.8 ±0.7 

8.36 ± 1.27 

Flat 

Flat 

Flat 

X-ray 

X-ray/SZE 

X-ray/SZE 


Note. — The results based on the combination of both weak and strong lensing are summarized 
(converted from quoted values assuming an NFW density profile if necessary). 

a Spherical: spherical prior on the intrinsic axis-ratios; Flat: flat prior on the intrinsic axis ratios; cos $: 
ACDM-like prior on the biased orientation of strong-lensing cluster halos (Corless et al. 2009). 
b External data sets used in combination with lensing for constraining the line-of-sight elongation. 
c The weak-lensing mass map of Umetsu & Broadhurst (2008) was used in the triaxial analyses by Oguri 
et al. (2005), Sereno & Umetsu (2011), Morandi et al. (2011), and Sereno et al. (2013). 
d NFW-equivalent of triaxial model parameters from Oguri et al. (2005). 



Figure 15. Marginalized constraints on the ellipsoidal NFW model param¬ 
eters (M 200 C 5 C 200 c) for A1689 compared to the c-M relations predicted 
for the full population of halos in ACDM cosmological simulations (Bhat- 
tacharya et al. 2013; Meneghetti et al. 2014; Diemer & Kravtsov 2015). The 
yellow shaded regions show the results from weak lensing combined with 
X-ray/SZE data. The red contours are from the full analysis of weak/strong- 
lensing and X-ray/SZE data. For each case, the contours show the 68.3% 
and 95.4% confidence levels in the c-M plane. The light blue areas show 
the la and 2cr ranges of intrinsic halo concentrations (with a 68% scatter of 
0.16 dex), respectively, as obtained by Diemer & Kravtsov (2015). All model 
predictions are evaluated at the cluster redshift z\ = 0.183. Overall, the in¬ 
ferred range of C 200 c is high but overlaps with the ~ 2 a tail of the predicted 
distribution for high-mass cluster halos. 

We have obtained evidence for a triaxial mass distribu¬ 
tion of A1689. The projected mass distribution derived from 
weak-lensing shear and magnification reveals a north-south 
elongation (ip e = 14.2° ± 8.4° east of north, see Table 5 
and Figure 1). We have determined the ellipticity of the pro¬ 
jected mass distribution to be e = 0.29 ± 0.07 (Table 5), 
which is typical for the population of collisionless CDM ha¬ 


los (Jing & Suto 2002) but slightly rounder than the stan¬ 
dard CDM prediction for the mean halo ellipticity, (e) ~ 0.4 
(Oguri et al. 2010). The matter ellipticity is detected at the 
4cr level from weak lensing alone, thanks to the greatly im¬ 
proved quality of Subaru data. Our free-form reconstruc¬ 
tion from HST strong lensing gives a consistent estimate of 
e = 0.27 ± 0.09. The ICM and matter distributions are co¬ 
aligned in projection (ip* = 12° ± 3°) but with different el- 
lipticities (ex = 0.15 ± 0.03), which is consistent with the 
geometric assumptions made (Section 6.2). 

When combined with X-ray/SZE observations, our lens¬ 
ing data favor a triaxial geometry of the matter distribution 
with minor-major axis ratio q a ~ 0.4 and major axis closely 
aligned with the line of sight (0 = 22° ± 10°, Table 7). 
These results are robust against the choice of priors and com¬ 
binations of lensing data sets. Despite that the intermediate- 
major axis ratio q /, is less constrained, the data prefer prolate 
( q a = qb) over oblate (qt = 1) configurations. A spherical 
configuration for A1689 is strongly ruled out. Overall, triaxial 
configurations fit the combined lensing and X-ray/SZE data 
much better than axially symmetric halos do (Sereno et al. 
2013). 

Our analysis shows that A1689 is elongated along the line 
of sight, as found by previous studies (Sereno et al. 2012, 
2013; Limousin et al. 2013). From the posterior samples, we 
find ey = 1.19 ± 0.37 (1.20 ± 0.34) and e[, CM = 1.22 ± 0.24 
(1.24 ± 0.25), as constrained by the combined weak/strong- 
lensing (weak lensing) and X-ray/SZE data sets. Such biased 
orientations are favored, although the intrinsic orientations are 
a priori assumed to be random. The a priori probability of a 
randomly oriented halo to have i? < 45° is ~ 29% (Sereno 
et al. 2013). The a posteriori probability of such a configu¬ 
ration is found to be 96% (99%) assuming a flat (7V-body- 
like) distribution of axis ratios. We emphasize that the use 
of X-ray plus SZE data is essential for obtaining data-driven 
constraints on the line-of-sight elongation. To break parame- 
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Table 10 

Ellipsoidal and spherically-enclosed mass estimates for A1689 


Overdensity a 

Ellipsoidal b 

Spherically enclosed 0 

A 

Ra M(< Ra ) 

r A 

±f sp h(< ri) 

500 

1.89 ±0.46 0.97 ±0.13 

1.08 ± 0.06 

0.88 ±0.13 

200 

2.79 ±0.69 1.24 ±0.16 

1.60 ±0.16 

1.15 ±0.16 


Note. — The overdensity radii are given in units of Mpc h 1 . The enclosed 
masses are in units of 10 15 Mq?i — L 


11 Mean interior overdensity with respect to the critical density p c for closure of 
the universe at z = 0.183. 

b Ellipsoidal overdensity radius Ra and total mass enclosed within Ra . 
c Spherical overdensity radius ta and spherically-enclosed total mass within r a ■ 

ter degeneracies in a lensing-only triaxial analysis, one would 
have to assume informative priors on the halo shape and ori¬ 
entation (Corless et al. 2009; Sereno & Umetsu 201 1). 

We find that the ICM is mildly triaxial with ~ 0.6 and 
c/[ CM ~ 0.7 (Table 8). The ratio of ICM to matter eccentrici¬ 
ties is e ICM /e = 0.87 ± 0.07 (Table 8), supporting the theo¬ 
retical assumption we have made that the shape of the gas dis¬ 
tribution is rounder than the underlying matter (Section 6.2). 

On the other hand, we find that the gas distribution is more 
elongated than the gravitational potential (e ICM /e A 0.7, Lee 
& Suto 2003), suggesting a deviation from HSE. These re¬ 
sults are again insensitive to the choice of the priors. The 
inferred values of (j/ CM and <^ CM are somewhat lower (more 
elongated) than, but consistent within errors with, the results 
of Sereno et al. (2012, 2013) based on the same X-ray data. 

The difference is mainly due to the improved, self-consistent 
SZE analysis. 


rfrxnc 

0.1 0.2 0.5 1 



Figure 16. Ratio of spherically-enclosed gas mass (M gas ) to total mass 
(Mtot) as a function of spherical radius r, derived from the full triaxial anal¬ 
ysis of weak/strong-lensing and X-ray/SZE data. The middle line tracks the 
median. The gray shaded regions represent the 68.3% and 95.4% quantiles of 
the distribution. Portions of these lines are dashed to indicate extrapolations 
to larger cluster radii. The horizontal bar shows the cosmic baryon fraction 
f b = determined by Planck Collaboration et al. (2015b). 


8.3. Gas Mass Fraction 


We compute the ratio of spherically-enclosed gas mass 
^ S ph,gas(< r) to total mass M sph , tQ t(< r) using the pos¬ 
terior samples of the ellipsoidal cluster model: 


/gas(< t) — 


-^sph,gas(^ f) 
M s p h,tot(< r) ’ 


(55) 


where M sp h(< r) denotes the total mass enclosed within a 
sphere of radius r, M sp h(< r) = J 4n dfl J Q r dr'r' 2 p(r') with 
did the solid angle. In Table 10, we list the values of ellip¬ 
soidal and spherical overdensity mass of the cluster evaluated 
at A = 200 and 500. 

The resulting / gas profile is shown in Figure 16 as a func¬ 
tion of integration radius r. The gas mass fraction within 
0.9Mpc ~ 1.2r- 2 5ooc is estimated as / gas (< 0.9Mpc) = 
O.lOO^QQ^g. When the gas mass measurements are ex¬ 
trapolated to r 5 ooc (Table 10), we find / gas (< r 500 c) = 
0.112+QQ20- When compared to the cosmic baryon frac¬ 
tion fb inferred from Planck Collaboration et al. (2015b), 
/gas (< r 500c )/fb = 0.7ll°;i2- These are consistent with typ¬ 
ical values observed for high-mass clusters (Allen et al. 2008; 
Umetsu et al. 2012; Okabe et al. 2014). 

Previous studies based on X-ray and lensing data found 
relatively low / gas values for A1689 using lensing total 
mass estimates, but assuming spherical symmetry: / gas (< 

0.25r 2O oc) = (0.0557 ± 0.0039)/i/ o 3/222 (Lemze et al. 
2008); / gas (< r 25 ooc) = 0.0552/“, / gas (< r 500 c) = 

0.0812/°;°^, and / gas (< r 20 0c) = 0.1053/°;°246 (Okabe 
et al. 2014, see also Kawaharada et al. (2010)). 

Umetsu et al. (2009) measured gas fractions for a sample 
of four high-mass clusters including A1689 from a joint anal¬ 
ysis of AMiBA SZE and Subaru weak-lensing observations, 
combined with published X-ray temperature measurements. 
Assuming spherical symmetry, they found for A1689 / gas (< 

r 2 500c) = 0.098/^26 and /gas (< r 500 c) = 0.115 ± 0.029, 
in excellent agreement with our results. Their gas fraction 
measurements are expected to be less sensitive to triaxiality 
because their / gas estimator depends on the ratio of the SZE 
and lensing signals, which are subject to similar projection 
effects albeit with somewhat different degrees of impact. 


R! r2ooc 



ICM ellipsoidal radius, /?[kpc/i '] 


Figure 17. Ratio of the thermal gas pressure (Pth) to the total equilibrium 
pressure (Ptot) in A1689 as a function of the ellipsoidal radius R measured 
along the major axis of the ICM halo. The middle line tracks the median. 
The gray shaded regions show the 68.3%, 99.4%, and 99.7% quantiles of 
the distribution, respectively. Portions of these lines are dashed to indicate 
extrapolations to larger cluster radii. 


22 Lemze et al. (2008) found r 200 c = 1-71 Mpc h 1 from their analysis. 
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8.4. Degree of Hydrostatic Equilibrium 

A quantitative assessment of the degree of equilibrium in 
the ICM is a critical issue for cluster cosmology based on 
hydrostatic mass estimates (e.g., Planck Collaboration et al. 
2014; Sereno & Ettori 2014). A significant advantage of 
our method is the ability to determine the intrinsic structure, 
shape, and orientation of the cluster system without a priori 
assuming HSE (Sereno et al. 2013). This allows us to com¬ 
pare the ICM properties directly to the gravitating mass cor¬ 
rected for projection effects, and thus to quantify the contri¬ 
bution of the thermal gas pressure P t h to the total equilib¬ 
rium pressure P tnt (Molnar et al. 2010; Kawaharada et al. 
2010). Here I\ ot is determined by the gravitational potential 
$ through VPtot = — Pg as V$ with p gas the gas mass den¬ 
sity. A consequence of the pressure equilibrium is the X-ray 
shape theorem (Buote & Canizares 1994), namely, that the 
gas in strict HSE is expected to follow iso-potential surfaces 
of the underlying matter distribution. For A1689, we find that 
the gas is more elongated than the gravitational potential (see 
Section 8.2), which points to a deviation from equilibrium. 

In Figure 17, we show the ratio of thermal to equilibrium 
gas pressure, Pth/Ptot, as a function of ellipsoidal radius It 
of the ICM distribution. For this aim, we have recomputed the 
posterior probability distributions for the cluster parameters, 
by imposing a sharp prior of e ICM /e = 0.7 (see Sereno et al. 
2013), corresponding to the assumption that the gas shape 
follows the gravitational potential. We find /%/Pot. ~ 0.6 
out to ~ 0.9 Mpc (~ 0.4r 2 oocX indicating a significant level 
(~ 40%) of non-thermal pressure support. The results here 
are consistent with Sereno et al. (2013), although our analysis 
favors a slightly higher level of non-thermal pressure support. 
We find no significant radial trend in the Pth/-Ptot ratio pro¬ 
file. 

Our results are in agreement with Molnar et al. (2010), who 
analyzed a simulated sample of massive regular clusters of 
(1 — 2) x 10 15 Mq!i - 1 having a smooth density profile, drawn 
from high-resolution cosmological simulations. Their simula¬ 
tions are therefore highly relevant to interpreting the observa¬ 
tions of A1689. They found a significant non-thermal contri¬ 
bution due to subsonic gas motions in the core region (20%- 
45%), a minimum contribution (5%-30%) at about 0.1r v j r 
(Lau et al. 2009), growing outward to about 30%-45% at the 
virial radius r v h (Nelson et al. 2014). 

Molnar et al. (2010) also tested the validity of HSE in 
A1689 using gravitational lensing (see Umetsu & Broadhurst 
2008; Kawaharada et al. 2010) and Chandra X-ray observa¬ 
tions under the assumption of spherical geometry, finding a 
non-thermal contribution of 40%. As discussed by Sereno 
et al. (2013), this however indicates that this test is highly sen¬ 
sitive to biases in the X-ray temperature measurements (Don¬ 
ahue et al. 2014). For the cluster, we find the Chandra tem¬ 
peratures are about 10% higher than the XMM-Newton results 
used here (Sereno et al. 2012), so that the thermal contribu¬ 
tion Pth/Ptot ~ 0.6 obtained by Molnar et al. (2010) could 
be correspondingly overestimated relative to our results based 
on the XMM-Newton temperatures. 

By combining Suzaku X-ray observations with the same 
lensing data as used in Molnar et al. (2010), Kawaharada 
et al. (2010) showed, assuming spherical symmetry, that the 
thermal gas pressure within r $ooc is at most 40%-60% of the 
equilibrium pressure and 30%^10% around the virial radius. 
Intriguingly, their Suzaku observations reveal anisotropic dis¬ 
tributions of gas temperature and entropy in cluster outskirts 


at rr,ooc, correlated with large-scale structure of galaxies 
surrounding the cluster. The outskirt regions in contact with 
low-density void environments have low gas temperatures and 
entropies, indicating that the outskirts of A1689 are in the pro¬ 
cess of being thermalized (Kawaharada et al. 2010). Their 
Suzaku temperature measurements are in agreement with the 
XMM-Newton results (Sereno et al. 2012). 

Morandi et al. (2011, see also Limousin et al. (2013)) ob¬ 
tained M 2 ooc = (1.81 ± 0.06) x 10 15 Mq/i _1 , c 2 oo c = 
5.71 ± 0.47, and q a ~ 0.5 for A1689 from a joint analysis 
of Chandra X-ray, weak-lensing, and strong-lensing data (see 
Table 9). The inferred level of triaxiality is similar to what we 
have found (Table 7), whereas the concentration is somewhat 
smaller and the mass is significantly higher than our results. 
They found that about 20 percent of the total ICM pressure is 
in non-thermal form, by assuming that Pth/Ptot is constant 
with radius and the gas shape follows the form expected for 
HSE. We note again that the Pth/Ptot results are also sensitive 
to calibration biases in the X-ray temperature measurements. 

The mass discrepancy between the present results and those 
by Morandi et al. (201 1) can be explained by the difference in 
their relative weights assigned to the weak- and strong-lensing 
data sets. As we have seen in Section 8.1.1, the HST strong- 
lensing data favor higher values of M(< r 20 o c ) (Table 7), 
although this represents a significant extrapolation beyond the 
radial range covered by the multiple images. Hence, if the 
parameter constraints are highly dominated by strong lensing, 
this could lead to an overestimate of M 2 oo c - 

8.5. Comparison with Planck data 

We compare the SZE measurements from the interferomet¬ 
ric data presented in Section 7.2 with a total power estimate 
based on the recent Planck data (Planck Collaboration et al. 
2015a). A1689 is detected by Planck with high significance 
(S/N > 15, Planck Collaboration et al. 2015d). We con¬ 
struct Planck SZE maps in two different ways with differ¬ 
ent assumptions, using the data in the 143 GHz, 217 GHz, 
and 353 GHz channels. The 217 GHz and 353 GHz bands 
are used primarily to remove the CMB and Galactic fore¬ 
grounds. The difference between the two maps accounts for 
different assumptions about the Galactic components: one 
is based on local estimates of the dust properties, and the 
other is on global properties. The resulting SZE maps are 
obtained at an effective resolution of 8' FWHM. The SZE sig¬ 
nal is integrated as a function of clustercentric radius. We 
obtain a direct estimate for the total Compton Y parameter 
of Ypi anc k = (3.8 ± 0.8) x 10“ 10 sr integrated out to a suf¬ 
ficiently large radius 13' (~ r 20 o c X beyond which the inte¬ 
grated SZE signal converges. Here the error is estimated from 
aperture photometry in the background regions 

This direct Planck measurement of the total SZE signal 
can be compared to the results inferred from the interfero¬ 
metric SZA observations (Section 7.2). Taking Ysza(< 6') 
(Table 6) as a lower limit on the total SZE flux, we find 
Ysza(< 6 ')/Ypianck = 1-45 ± 0.44. Hence, the results from 
two independent SZE instruments operating at different angu¬ 
lar scales are compatible with each other at la. The relatively 
low Y value derived from the Planck data could be understood 
in light of the low gas temperature and entropy at £ r-.ooc ob¬ 
served by the Suzaku X-ray satellite (Section 8.4). The Suzaku 
X-ray observations are in agreement with the thermal pressure 
profile of A1689 obtained from Planck data out to ~ 2r5oo c 
(Y. Mochizuki et al. 2014, submitted to ApJ). 

When compared to Planck’s hydrostatic mass estimate. 
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M 500 C = (8.77 ± 0.34) x 1{) [A MqKjq , our lensing mass 
measurements (Table 10) give a spherical mass ratio of 
Mpianck/M gl = 0.70±0.15 and 0.58±0.10 with and with¬ 
out corrections for lensing projection effects, respectively. 

9. SUMMARY 

We have carried out a 3D multi-probe analysis of the rich 
cluster A1689, one of the most powerful known lenses on the 
sky (0Ein = 47.0" ± 1.2" at z s = 2, Table 1), by combining 
improved weak-lensing data from new wide-held BVRci'z' 
Subaru/Suprime-Cam observations (Sections 3 and 4) with 
complementary strong-lensing (Section 5), X-ray and SZE 
(Section 7.2) data sets. 

We have generalized the ID weak-lensing inversion method 
of Umetsu et al. (2011b) to a 2D description of the mass dis¬ 
tribution without assuming particular functional forms (Sec¬ 
tion 2). This free-form method combines the spatial shear 
pattern with azimuthally averaged magnification information, 
the combination of which breaks the mass-sheet degeneracy. 

We have reconstructed the projected matter distribution 
from a joint weak-lensing analysis of 2D shear and az¬ 
imuthally integrated magnification constraints (Section 4). 
The resulting mass distribution reveals elongation with an axis 
ratio of q± ~ 0.7 in projection (Figures 1 and 8), aligned well 
with the distributions of cluster galaxies and ICM (see Kawa- 
harada et al. 2010). When assuming a spherical NFW halo, 
our full weak-lensing analysis yields a projected halo concen¬ 
tration of clooc = 8-9 ± 1.1 (c^P ~ 11), which is consis¬ 
tent with and improved from earlier weak-lensing work based 
on Subaru Vi' imaging (Umetsu & Broadhurst 2008; Umetsu 
etal. 2011b). 

We obtain excellent consistency between weak and strong 
lensing in the region where these independent data overlap, 
200 kpc (Figures 6 and 10). We also find an improved 
agreement between weak and strong lensing in terms of con¬ 
straints on projected NFW parameters (Figure 11) relative to 
previous work (Sereno & Umetsu 201 1). This is largely due to 
improved techniques for strong-lensing reconstruction and to 
careful regularization of the covariance matrix (Section 7.1.2). 

In a parametric triaxial framework, we have determined the 
intrinsic stmcture, shape, and orientation of the matter and 
gas distributions of the cluster, by combining weak/strong 
lensing with X-ray/SZE data under minimal geometric as¬ 
sumptions (Section 7). We have shown that the data fa¬ 
vor a triaxial geometry with minor-major axis ratio q a = 
0.39 ± 0.15 and major axis closely aligned with the line 
of sight (i9 = 22° ± 10°). A spherical configuration for 
A1689 has been strongly ruled out. We obtain a halo mass 
M 200 C = (1-24 ± 0.16) x 10and a halo concentra¬ 
tion C 200 C = 8.36 ± 1.27, which is higher than typical con¬ 
centrations found for high-mass clusters (3 .1) C 200 C ~ 6; e.g., 
Okabe et al. 2013; Umetsu et al. 2014; Merten et al. 2014), but 
overlaps well with the £ ler tail of the predicted distribution 
(Figure 15; Bhattacharya et al. 2013; Meneghetti et al. 2014; 
Diemer & Kravtsov 2015). 

We find that the ICM is mildly triaxial with q]' : M = 
0.60 ± 0.14 and ql GM = 0.70 ± 0.16 (Table 8). The gas 
distribution is rounder than the underlying matter, e ICM /e = 
0.87 ± 0.07, but more elongated than the gravitational po¬ 
tential (e ICM /e >, 0.7), suggesting a deviation from equilib¬ 
rium. The gas mass fraction enclosed within a sphere of radius 
r = 0.9Mpc ~ 1.2r25ooc is found to be / gas = 10.0/f'g%. 
When the gas mass measurements are extrapolated to r$ ooc, 


/gas(< T 5 ooc) = H.2/^%. When compared to the cos¬ 
mic baryon fraction f b (Planck Collaboration et al. 2015b), 
we find / gas (< r 500 c)/fb = 0.71(Figure 16). These are 
consistent with typical values observed for high-mass clus¬ 
ters. The thermal gas pressure contributes to ~ 60% of the 
total pressure out to ~ 0.9 Mpc (Figure 17), indicating a sig¬ 
nificant level of non-thermal pressure support. The results are, 
however, sensitive to calibration biases in the X-ray temper¬ 
ature measurements (Donahue et al. 2014). When compared 
to Planck ’s hydrostatic mass estimate, our lensing mass mea¬ 
surements yield a spherical mass ratio of Mpi anc k/ Mql = 
0.70 ± 0.15 and 0.58 ± 0.10 with and without corrections for 
lensing projection effects, respectively. 

Extending this work to larger samples of clusters will en¬ 
able us to recover intrinsic distributions of cluster structural 
properties (e.g., M 20 oc, c 2 ooc) and axis ratios (q a , q b ), for a di¬ 
rect statistical comparison with the standard ACDM paradigm 
and for a wider examination of alternative DM scenarios (e.g., 
Schive et al. 2014). The CFASH survey (Postman et al. 2012) 
provides such ideal multiwavelength data sets of high quality 
(Donahue et al. 2014; Umetsu et al. 2014; Zitrin et al. 2014; 
Czakon et al. 2014; Rosati et al. 2014), for a sizable sample 
of 25 high-mass clusters. 
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APPENDIX 


A. NONLINEAR EFFECT ON THE SOURCE-AVERAGED LENSING FIELDS 
A.l. Reduced Gravitational Shear 

The reduced shear, g = 7/(1 — re), is nonlinear with re, so that the averaging operator with respect to the source redshift 
acts nonlinearly on re. In general, a spread of the source redshift distribution, in combination with the single source-plane 
approximation, may lead to an overestimation of the gravitational shear in the nonlinear regime. 

Let us expand the reduced shear g = g(z) with respect to k(z) = W(z)n oo and 7 (z) = W(z) r ) 00 as 

OO 

g = 7/(1 - re) = W loo ( 1 - ILreoo )- 1 = W loo (W Koo ) k . (Al) 

k—0 

The reduced shear averaged over the source redshift distribution is expressed as 


(g) = 

k =0 


(A2) 


where the angular brackets represent an ensemble average over the redshift distribution of background sources. In the weak- 
lensing limit where k x <C 1, (g) ~ (W) g r y 00 = ( 7 ). The next order of approximation is 


(d)g ~ Too ((W)g + fW )g K 00) 


(W ) gl oo 

l-Koc(W 2 ) g /(W) g - 


(A3) 


Seitz & Schneider (1997) showed that Equation (A3) yields an excellent approximation in the mildly-nonlinear regime with 
re 00 ~ 0.6. Defining fw, g = (W 2 ) g / ( W) 2 g , we have the following expression for the source-averaged reduced shear valid in the 
mildly-nonlinear regime: 


(5) 


(7) 

1 - fw,g(K) ’ 


(A4) 


with (re) = (W) g K 00 . For a lens at relatively low redshift, fW 2 ) g « (W) 2 and fw, g ~ 1, leading to the single source-plane 
approximation: (g) rj ( 7 )/(l — (re)). The level of bias introduced by this approximation is A g/g ss (fw,g ~ l)(re). In typical 
ground-based deep observations of z\ x 0.5 clusters, A fw — fw — 1 is found to be of the order of several percent (Umetsu et al. 
2014), so that the relative error is negligibly small in the mildly-nonlinear regime. 


A.2. Magnification Bias 

Let us consider a maximally-depleted sample of background sources with a = —d log N t ,(> F)/dF = 0, for which the effect 
of magnification bias is purely geometric, b g = /z -1 , and insensitive to the intrinsic source luminosity function. In the nonlinear 
subcritical regime, the source-averaged magnification bias is expressed as (Umetsu 2013; Umetsu et al. 2014) 

(/z- 1 ) = (1 - (re )) 2 - |( 7 )| 2 + (fw,, ~ 1 ) ((«) 2 - ( 7 ) 2 ) « (1 - («)) 2 - l( 7 )| 2 , 


(A5) 
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where fw.n = (W 2 )n/ (W - )^ is of the order of unity, (k) = (W) tl K 00 , and ( 7 ) = (W)^ 7 oo- Hence, the error associated with 
the single source-plane approximation is (A/ 1 -1 ) = (fw.n ~ 1)((k ) 2 — (t) 2 ) = — (l) 2 ), which is much smaller 

than unity for background populations of our concern (A fw.n ~ O(10 -2 )) in the mildly-nonlinear subcritical regime where 
(k) ~ | ( 7 ) | ~ O(10 _1 ). It is therefore reasonable to use the single source-plane approximation for calculating the magnification 
bias of depleted source populations with aCl. 


B. DISCRETIZED EXPRESSIONS FOR CLUSTER LENSING PROFILES 

First, we derive a discrete expression for the mean interior convergence (< 8 ) as a function of clustercentric radius 8 using 
the azimuthally averaged convergence (O'). In the continuous limit, the mean convergence k x (< 8 ) interior to radius 8 can 
be expressed in terms of ( 8 ) as 

9 r e 

Koo {<0) =J dln8'8 ,2 Koo (8'). (Bl) 

For a given set of (iVbin + 1) concentric radii 8 i (i = 1,Abi n + 1), defining 7Vb; n radial bands in the range O r nin = 81 < 8 < 
0 jv bin _l_i = 0 max , a discretized estimator for k 00 (< 8 ) can be written in the following way: 


ttoo 


0 z 


2 

Noo(^ ^min) T 


2 i_1 

sEAln^Voo^i), 

i j =1 


(B2) 


with Aln 0 j = ( 0 j + 1 
we have 


8 i)/ 8 i and 8 i the area-weighted center of the ith annular bin defined by [0*, 0,; + i]. In the continuous limit. 


8i = 2 



dd' 8 , 2 /{ 8 2 +1 


2 8 2 + 8 2 +1 +8A+i 

3 d.i + 8 i+ i 


Next, we derive discretized expressions for the tangential reduced shear g+{ 8 ) and the inverse magnification 7 1 (0) in terms 
of the binned convergence n 00 ( 8 i ), using the following relations: 


\ _ (W)g [ K oo(< $i) ^oo(^z)] 

9+[ i} ~ i-fw, g (w) gKoo (di) ’ 

d~ 1 (0 l ) = [ 1 - (W)„ Koo ( 9 z )] 2 - (W)l [ Koo (< 9i) - Koo^i )] 2 , 


(B3) 

(B4) 


where both the quantities depend on the mean convergence interior to the radius 0,. (< Si). By assuming a constant density 

in each radial band, we find the following expression for «oo(< 0 ?): 


Koo(< Oi) = ^ (di/di) 2 Koo(< 9i) + (9 i+1 /9i) 2 Koo(< 8i+ 1 ) 


(B5) 


where Koo(< 8 t ) and Koo(< 0i+i) can be computed using Equation (B2). 

Accordingly, all relevant cluster lensing observables, g+(8) and n^(8), can be uniquely specified by the binned convergence 

profile {ftoo,min: 1 with tv 00 ,min = 0min) and hioa.i = K'oo(^z)* 


C. TWO-DIMENSIONAL TO ONE-DIMENSIONAL PROJECTION 

To make a direct comparison between the results from ID and 2D weak-lensing analyses, we construct a projected mass profile 
£(0) from an optimally weighted radial projection of the E(0) field as (Morandi et al. 201 1) 


E(i) = 




A t C7J'E 


( 2 )^( 2 ) 


(Cl) 


where £( 2 ) = {£(0m)}m =1 i s a pixelated mass map, C( 2 ) is the pixel-pixel covariance matrix of 5T 2 ), 2^) is a vector of 
radially binned £ values, and A is a mapping matrix whose elements A,„, represent the fraction of the area of the mth pixel lying 
within the ith clustercentric radial bin (Section 2.3.2). The covariance matrix for E (1 ) is given by 


C'd) = 




-1 


(C2) 











